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ABSTRACT 


One-dlmenslonal quasi-steady theory Is used to develop an engineer- 
ing method to determine the time-dependent pressure In vented launch- 
vehicle compartments during the first few minutes of the ascending 
flight. 

This method, by nature an Iteration procedure. Is Intended to pro- 
vide the structural design engineer with the reduced loads on the com- 
partment walls resulting from the venting process. The basic program Is 
set up for the CDC-3200 digital computer which can handle presently only 
up to N ® 5 compartments where the Inner compartments must have only 
one connecting orifice and the last compartment can have up to NV = 5 
orifices venting Into the atmosphere. Furthermore, the compartments 
have to be placed In series. Though the compartment and orifice number 
can be raised Indefinitely, It Is advisable to restrict the number to 
as few as possible to keep the computation time low. 

The basic program has been extended to offer combinations of com- 
partment and connecting orifices. Compartment leaks and their accompany- 
ing coefficients, as well as venting through a duct of varying cross 
section, have been Included. The effect of heating or cooling the duct 
flow can also be computed. 
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DEFINITION OF SYMBOLS 


Symbol 

Definition 

Dimension 

A 

orifice area 


a 

mean radius of earth 


C 

radial clearance 

m 

c 

velocity of sound 

m sec'^^ 

H 

geopotentlal altitude 

m 

L 

flow length 

m 


hydraulic radius, diameter 

m 

V 

volume 


X 

length of tube 

m 

Z 

geodetic altitude 

m 

®o 

gravity 

m sec“^ 

M 

Mach number 

- 

t 

time 

sec 

V 

velocity of gas 

m sec"^ 


specific heat at constant pressure 

m Kg Kg"^ 

Cy 

specific heat at constant volume 

m Kg Kg*^ °K 

f 

friction factor 

- 

h 

enthalpy 

m Kg 

K. ^0 

loss coefficient, contraction coefficient 

L' 

temperature gradient 

°K Km-^ 

m 

mass 

Kg sec^ m"^ 

m 

mass flow rate 

Kg sec m"^ 


iv 


DEFINITION OF SYMBOLS 


Symbol 

P.P 

Q 


R 

T 

w. 

7 

P 

P 


0 


Definition 

molecular weight 

pressure, mean pressure 

heat per unit mass of gas entering 
the volume 

Reynolds number 

gas constant 

temperature 

external work 

specific heat ratio 

density 

absolute viscosity 
shearing stress 

loss coefficient for pipes other than 
friction 

area ratio Aq/A^ 


Dimension 
Kmol Kg"^ 

Kg m“2 

Kg m® Kg"^ sec"^ 


m Kg Kg"^ °K"^ 
®K 

Kg m 

Kg sec^ m"'* 

Kg m”^ sec“^ 

Kg m-2 


Subscripts 

1,2, i designates sections under investigation 

t total flow properties, total pressure, total density 


Throughout this report the technical system of units is used; therefore 
the unit Kg is the Kg force = Kp (Kilopond). 
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COMPARTMENT VENTING AND PIPE FLOW WITH HEAT ADDITION 

SUMMARY 


One-dimensional quasi-steady theory is used to develop an engineer- 
ing method to determine the time-dependent pressure in vented launch- 
vehicle compartments during the first few minutes of the ascending 
flight. 

This method, by nature an iteration procedure, is intended to pro- 
vide the structural design engineer with the reduced loads on the com- 
partment walls resulting from the venting process. The basic program is 
set up for the CDC-3200 digital computer which can handle presently only, 
to N = 5 compartments where the inner compartments must have only one con- 
necting orifice and the last compartment can have up to NV = 5 orifices 
venting into the atmosphere. Furthermore, the compartments have to be 
placed in series. Though the compartment and orifice number can be raised 
indefinitely, it is advisable to restrict the number to as few as possible 
to keep the computation time low. 

The basic program has been extended to offer combinations of com- 
partment and connecting orifices. Compartment leaks and their accompany- 
ing coefficients, as well as venting through a duct of varying cross 
section, have been included. The effect of heating or cooling the duct 
flow can also be computed. 


I. INTRODUCTION 


Any internal launch-vehicle space bounded by rigid walls is defined 
as a compartment. These compartments must be depressurized or vented 
during the ascending portion of the flight trajectory when the surround- 
ing atmospheric pressure decreases rather rapidly. Depressurization is 
generally obtained by orifices cut into the rigid wall, or openings 
leading into another compartment which is vented to the surrounding 
atmosphere. The size and shape of the hole, as well as the pressure 
differential from one chamber to another or to the atmosphere, deter- 
mine the mass flux out of the compartment. Because of the mass transfer 
to the next lower compartment, the density in the first compartment will 
decrease and so will the internal pressure. Any inflowing mass from 
upstream compartments, however, will increase the pressure again. This 


indicates the dependence of the gas properties of the compartment on 
downstream and upstream conditions. However, should choking occur In 
one of the orifices or In the vent duct, the line of dependence Is then 
restricted to the downstream or upstream portion of the flow only. For 
the Isotropical orifice flow, empirical correction factors are used 
which are a combination of contraction coefficients (vena contracta), 
friction coefficients, and, for the outside vents, the Influence of 
the outside tangential flow. 

Occasionally compartments are vented through ducts which might con- 
sist of pipes with varying diameters. Adiabatic flow theory Is then 
applied with suitable pipe friction factors to calculate flow properties 
In the duct. Heat addition to the flow Is neglected when the duct Is 
not too long, since the time range involved is relatively short. How- 
ever, for long pipes heat addition to the flow can have a definite effect, 
especially when the temperature difference between ambient and pipe flow 
is large. 

The equations for the atmospheric data are taken from the U. S. 
Standard Atmosphere, 1962 [9] with a simplified calculation of the geo- 
potentlal altitude. These data can be corrected for a specific launch 
location if a table of correction factors for the pressure Is added as 
a subroutine to the program. 


II. THE ANALYTICAL APPROACH 


A. Simple Orifice Flow 

The venting problem, in reality a time-dependent unsteady 
problem, is simplified to a quasi-steady one -dimensional problem; all 
properties are therefore uniform over each cross section. Further 
assumptions are that the heat transfer between the surrounding com- 
partment walls and the gas within is negligible, and the gas remains 
perfect throughout the entire venting process. 


The approach consists of applying the mass and energy balance 
equation and the equation of state. We begin with the mass conservation 
law, expressed in general as 




a(m) 


cv 


St 




( 1 ) 


which states that the rate of accumulation of mass within the control 
volume (c.v.) is equal to the excess of the incoming rate of flow 



over the outgoing rate of flow 

1 *o- 

At any instance of flow, 



m cv 


where dV is an element of the control volume, p is the local mass density 
of the control volume, and the integral is to be taken over the entire 
control volume. Furthermore, 





CV 


bt 


dV 


and 







I 


pVn dAo. 


(3) 


(4) 


Equation (1) can now be written as 


/ It - f P^n ^o. (5) 

cv 


3 



where Vn is the velocity component normal to dA. The control volume is 
now supposed to be the compartment, and with the assumption of one- 
dimensional flow we obtain 


3f ’ ■ Z *1 ■ I *»• 

' n m 


For steady flow, dp/dt = 0, and the incoming mass is equal to the out- 
going mass. In equation (6), n and m designate the openings in the com- 
partment through which gas can flow into and out of the chamber, 
respectively. For the venting problem, we allow the density p to change 
with time according to the difference of mass flux rate. 

The energy equation of steady flow relates the external work 
effect and the external heat exchange to the increase in the flux 
enthalpy, kinetic energy, and potential energy passing through the con- 
trol surface; then 


m(dQ - dWjj) = m(h + dh) - mh + m + d - m ^ + m(z + dz) - mz, 

(7) 

where dQ is the net heat added to the stream from sources external to 
the main stream per unit mass of gas entering the control surface. 
Likewise, dWj^ is the external work delivered to the outside body per 
unit mass of gas entering the control boundary. For our problem, dW^ “ 
and z, the height of the stream centerline to datum line, is negligible 
for gas flow. Then equation (7) becomes 


mdQ = mdh + md(v^/2). 


( 8 ) 


The enthalpy h is a function of the temperature only; therefore. 


dh = CpdT, 


(9) 


4 



and equation (8) becomes 


mdQ = mCpdT + md(v®/2). (10) 

According to our assumption that no heat is added from the outside, 
dQ = 0, and equation (10) becomes 

dh = -d(v2/2). (11) 

The mass flow rate is now given by the continuity equation as 

m = pvA. (12) 

With the equation of state 

P = pRT, (13) 

and 
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the velocity can be obtained from the energy equation for adiabatic flow, 


+ 2c T = const. 

P 

By assuming that the velocity in the upstream compartment is negligibly 
small, Vi = 0, we obtain finally for the velocity with which mass is 
flowing out of the compartment (Figure 1) 


V 


ex 



r P n 
Pi ^ ex 

Pi P 

^ex-J 


(14) 


5 



I 


with the relation between pressure and density for an Isentroplc process 
of a perfect gas between two stations, 


^ ( 15 ) 


substitution of equations (13) through (15) Into (12) yields the mass > 

out of the compartment: 


m = AK ^/2PlPl 




» 

(16) 


The discharge from a compartment orifice Is considered to be an Isen- 
tropic process, since we used isentroplc relations for establishing the 
mass flow rate. Any losses of mass flow rate due to total pressure and 
contraction are represented by the loss coefficient K < 1. The pres- 
sures Pex snd Pi, as well as the density p^, have to be representative 
mean values for a time-dependent process. 

If the flow In the orifice becomes sonic, the mass flow Is 
Independent of the pressure difference across the orifice: 

1 

A - AkJ2Pipi /t?T- 


Because of the outflowing and inflowing mass, the density, pressure, and 
other properties within the compartment are changed. To obtain the 
necessary relation from equation (6), we approximate the differential 
quotient dp/dt by a difference quotient: 


dp ^ Ap ^ P(t + At) - P(t) p(ti) - p(to) 
dt At At ti - tjj 


(18) 
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Substituting equation (IS) into equation (6) results in the following 
expression for the density: 


At 


P(ti) “ p(to) 


1 + 




p(t«) V 




(19) 


and according to the gas relation (equation (15)) the pressure at time 
ti becomes 


P(ti) = P(to) 



( 20 ) 


The pressure, density, and mass flow terms in equations (20) and (16) 
have to agree with each other within the compartment during a particular 
time step. At * (t^ - tp) . Since equation (16), the mass flow rate equa- 
tion, uses mean values of the pressures, an Iteration procedure is neces- 
sary so that Pi represents a mean value of P(ti) and P(t{j) in the particular 
time interval At. 

To circumvent at least one Iteration procedure for finding the 
mean pressure value of a single compartment, a series approximation of 
equation (16) and (20) may be used provided that the pressure change for 
a small time step can be assumed to be linear. 

If we set 


Pi(ti) + Pi(to) 
Pi = Pi o 


and if 


P 


ex 


P 


1 


7 



(where Pgj^ is also a mean value of the external pressure at that particular 
time step) does not differ very much from unity, a series approximation 
can be made with 


X 


h ^ 

Pi(toV 


« 1 . 


(21) 


Departing from equation (16) with the above notation, we obtain for n 
vent holes 



With the isentropic gas relation, we can introduce the pressure of the 
compartment at time (tp), P(to), which is assumed to be known since it 
is the pressure at t^ of the previous time step. 


P P 

ex ex 


Pi(to) 


Pi Pl(to) Pi 


and the square root term is 


n/ 2PiPi = N/2Pi(t^) pi(to) 



y+1 


(23) 


(24) 


With the abbreviations 


*^o n/ y - 1 *^2Pi(to) Pi(tQ) 



1 


VPi(to) 





(25) 
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and 



( 26 ) 


we obtain for the mass flow rate of equation (22) finally 





(27) 


and with equation (21), equation (27) becomes 


_2 



The first term on the right-hand side of equation (28) can be developed 
into a series 


1 - X 


7 


■ 1 - 




Inserting (29) into (28) yields 



Cl - X - x2 - ~ ^ X3 + O(X^). 

( 30 ) 


9 



similarly, we proceed with equation (20), where the outflowing mass Aq is 
brought to the left side. With the stipulation already mentioned above, 
i.e. , 


_ P(ti) + PCtg) 

p . 


we obtain 


I 


A- 



+ 


Pi(to)V 

At 


1 



(31) 


With the substitution of 


1 



1 + — ^ ^ 
Pl(to)V 


(32) 


and equation (21) inserted into equation (31) yields 


— 2 


At Cq 


Z 




= 1 - 2C; 



1/7 

- 

2(1 - X) - 1 

+ c 

2 

2(1 - X) - 1 


2/7 


(33) 


The inflowing mass 



in equations (31) and (32) is assumed to be known, since it represents 
the mass-flow of the upstream compartment for which the procedure described 
here has already been conducted. For n(n > 1) compartment, however, a 


10 


second iteration is necessary, since the outflowing mass, mo, influences 
the inflowing mass, m£, unless the upstream orifice was choked. Develop- 
ment into a series leads finally to 


2 


At C 2 


I '"c 


Pi(to) V 


(1 - Cs )2 + 


4cpX 


(1 - ca) 


+ i£ 2 Xf 


(7 - 1) + C2(2 - 7) 


(7-1) 8 c 2 X=^ 

+ 373 


(27 - 1 ) + 4ca(2 - 7 ) 


+ 0(X^). (34) 


Equations (34) and (30) have to be set equal in iho to obtain X. Before 
we go to this step, however, some useful abbreviations are introduced. 
We set; 


0 ^ = (1 - Ca) 2 ; 

= ^ca(l - Cs)/ 7 ; 

0)2 - 4c 2[(7 - 1) + ca(2 - 7)]/7®; 

0(3 = (7 - 1) 8 c 2[(27 - 1 ) + 4ca(2 - 7)3/37^; 

0:4 = [AtC 2 /(Pi(to)V)]'" 


and 


3' * 
3 


Pq/ P4; 

^0 “ 

Pi/ P4 ; 

a'l = o:i/a4 

Pa/ P4 » 

a' = a2/Q!4 

P3/ P4; 

= a3/a4 


Po “ 

(1 

“ Cl) 

Pi - 

(7 

- l)/7 

Pa ■ 

(7 

- 1 )/( 272 ) 

P3 - 

(7 

- 1)(7 + 1)/(67^) 

P4 “ 

(1/ 



A • . B'i + Oil 

1 PL + aL 


p' - a,' 
^o o 

PL + 


a' s 0^3 + p3 

3 p' + a’ 


(35) 
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and the solution for X yields 


X = - ^ + n/a; - X3 + (a;^/2)2 , (36) 

where A 3 can be left out and X can be obtained without Iteration or can 
be used as a correction to X, If retaining of only the quadratic terms 
Is not enough to obtain an adequate accuracy. However, only the positive 
sign of the square root is used for the mass flowing out of the compart- 
ment. The mean pressure of the compartment is finally 




Pi “ Pi(to)[l - X]. 


(37) 


For choked flow, equation (17) is used, which becomes now, for n orifices 
of the compartment. 


2±1 


_Z " Gife) Z C + 0 


7-1 


n 


With the abbreviation 


7 + 1 


2Pi(tQ) Pi(to) 


]• 


(38) 




2Pi(to) 



I 




(39) 


equation (38) becomes 



2~±j. X + X^ + X^ 

7 272 673 


+ O(X^) . 


( 40 ) 
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Substituting 



9 ^ 

n 

o 

CCL 

Po " Po/P4 


Pi - (7 + i)/7; 

Pi - P1/P4 


P2 “ -(7 + 1)/(27^); 

pL - P 2/P4 

% 

Ps * -(7 + 1)(7 - 1)/(67^); 

P3 ■ P3/ P4 


P4 “ 1/c*; 



we obtain the same expression as equation (36). 


B. The Adiabatic Flow in Constant Area Ducts with Friction 


The assumption is again that the gas is perfect. The rate of 
change of the gas properties depends now on the amount of friction, so 
that the momentum equation must be introduced. 

The perfect gas relation was given by equation (13); taking 
logarithmic differentials, we obtain 


^ le 4 . 

P " p T • 


(41) 


The definition of the Mach number, which was not previously introduced 
is, for a perfect gas. 


= v^/yRT 


(42) 


\ 



dv‘ 


or 


T • 


(43) 


The energy equation, given by (10) with dQ = 0 becomes by dividing 
through by CpT and using the definition for the Mach number 


T 




0 . 


( 44 ) 


The continuity equation (12) becomes in terms of logarithmic differen- 
tials with the mass flow as a constant 


^ + i dvf . „ 

p 2 v*' 


The momentum equation is written as 


-AdP - tjj dAj^ = Adv, 


(45) 


(46) 


where Tjj is the shearing stress exerted on the stream by the walls, and 
dAx is the wetted wall area over which acts. Introducing a friction 
coefficient f(x), we obtain 

2x 

f = -I (47) 

pv=^ 


and a hydraulic diameter 


_ 4A _ fs 

" dA /dX ■ 2 * 

X 


(48) 


For a circular pipe the hydraulic diameter is thus equal to the diameter 
of the pipe. Inserting these expressions into equation (46), we obtain 
the momentum equation which, divided through by P, can be written, with 


P ^ 2 Djj 2 ^ 
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By eliminating dP/P and dv^/v®, we obtain finally, for a pipe of length 
L between sections 0 and L 


dM^. (49) 


With a mean friction factor f between stations 0 and L as 




Equation (49) when integrated between 0 and L or M§ and M| becomes 



Known data are the friction coefficient f for the pipe of length L and 
the hydraulic diameter D^. The inlet Mach number m| can be computed 
with the aid of the mass flow from the upper compartment to the pipe 
and the pressure at the pipe inlet, depending on the type of approach: 
isentropic for a converging pipe inlet section or adiabatic for a sudden 
constriction. The pressure, not known a priori, can only be obtained by 
an iteration scheme, assuming first a Mach number, Mq. Then with equa- 
tion (51) the Mach number increase for subsonic flow can be calculated 
between 0 and 1. With the relation 



15 



the respective pressure drop due to friction can be calculated. If the 
Mach number at station 1 is subsonic, then the exit pressure Px must be 
equal to Pex> the atmospheric pressure outside the pipe. If this is not 
the case, we must repeat the calculation with another Pg, which in turn 
will alter the mass flow from the upstream compartment. 

If the exit of the pipe is choked, Mi “ 1 , then the Mach number 
Ml is fixed and one has to calculate on an iterative basis the inlet 
Mach number. Mg, and the corresponding pressure, Pg. The pressure and 
density at station 0 are then given by the relations 


P 


o 


p* 



y + 1 



(52) 


and 


V 

o 


V* 




y + 1 



(52a) 


where the asterisks designate critical values for which the Mach number 
is M = 1. 


Equation (51) cannot be solved explicitly. For subsonic flow, 
the following iteration scheme has proven to be stable. 


Ml = 


M 


o 



(53) 


4 




where ^ is a loss coefficient of obstacles in the pipe such as bends, 
constrictions, etc. It can be a function of the Mach number or the 
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geometry of the pipe. In the equation Mi appears in the logarithmic 
term of the denominator. The length of the pipe for choked flow is 
given by 


(1 + 7 ) 

o 


C. Flow in Ducts with Heating or Cooling and Friction 

The determination of the flow parameters again requires the 
application of the three conservation equations and the perfect gas 
equation. We retain now dQ in equation (10): 


dQ = c dT + d(v^/2) = c dT^., 


and with the definition of the Mach number of a perfect gas equation 
(42), we can find the following relations [8]: 

For the Mach number increase along the pipe, 


iHF Tt(x) 1 - (Dh/L) ^ ^ 


For the velocity increase, 


, 1 + ^ ~ M^ dT^ 

dv 2 t 


2 t yM^ . ? d(x/D 

1 - m 2 T^(x) 2(1 - m 2 ) ^ (Djj/D* 


For the speed of sound. 


;2) ( 1 + m 2 ) dT. 


(1 - 7M^) [} 

2(1 - M^) 


T^(x) 4(1 - m2) 


^4^4f 4^. 


(Djj/L) 



For the temperature, 


^ (1 - (l + ^ M^) dij 

T “ 1 - T^(x) 


y(y - 1) - d(x/L) 

2(1 - (Djj/L) 


(59) 


For the density decrease along the pipe, 


P 


1 + dT^ 

7S‘ 




T^(x) 


4f d(x/L)., 

2(1 - M^) (Dj^L) 


(60) 


For the pressure decrease, 


7(m2) (l + dT^ 

T “ 1 - Tj. (x) 



( 61 ) 


The total pressure loss Is given by 


= _ 2^ 4 f d (x/L) . 

2 T^(x) 2 Td^ 


(62) 


Equations (56) to (62) have to be integrated numerically. Q is the heat 
transferred to the flow within the tube from external sources, per unit 
mass of stream. 

Equation (56) shows that heat addition, expressed as total 
temperature increase, causes the Mach number to increase above that 
caused by the friction force. The choking length of the tube is there- 
fore shorter. We notice that equation (56) is a singular ordinary dif- 
ferential equation of the first order. In the domain 0 S < 1, the 
right-hand side of equation (56) is continuous. A singular point is 
encountered at = 1. 
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We solve the equation by applying the method of successive 
approximation, which Is also known as Picard's method. If equation (56) 
Is written In the general form as 


dM^ 

dx 


f (M^,x), 


then the solution has the form of 


X 

M^(x) = M^(Xq) + J d|. 


This relation Is, In reality, an Integral equation. Involving the dependent 
variable under the Integral sign. It can be solved on an Iterative basis 
as follows: 


X 

m2(x) -m2(Xo) + J f(|,M§) dg 


X 

M|(x) =m2(x^) + J f(|,M2) d| 


X 

M^(x) = m2(Xq) + J f(g,M^.^) dg. 


As n Increases Indefinitely, the sequence of functions M^(x) tends to a 
limit which is a continuous function of x, and the limit- function satis- 
fies the differential equation. On the other hand, the Lipschitz con- 
dition has to be satisfied. If (x,M^) are two points within the domain 
of the same abscissa, then 
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f(x,M|) - f(x,M2)| < k1m| - Mfl, 


where K is a constant. In the neighborhood of M = 1, the Lipschitz con- 
dition no longer holds. Rearranging equation (56) yields 


dx 

dM^ 



1 - 


(1 + 


dT. 


d(j;/L) T^(x) 


+ 7M' 


2 4f h 


or in general form 


dx 

dM‘= 


= g(x,M=). 


The function g(M^,x) is now singular at = 0 and dx/dM^ = 0 for = 1. 
The method of successive approximations can again be applied. The solu- 
tion now has the form 




x(M^) = 


x(M§) 


g(x,M^) dM^ 


Mf 


Instead of the Lipschitz condition, a less stringent condition can be 
applied. For instance, a step sensitivity condition such as 


i dM" 
^ dx 




e^|xi - X, 


for the first case and 


dx 

dM“= 


s £2; 


xi - XqI g ■ ^o 
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If and 62 are of the order €i;e2 ^ .002, the analytic solution of 
equation (56) with dT^. = 0 was almost exactly reproduced. The two solu- 
tions are joined where the step sensitivity condition was satisfied for 
both regions. 

Once the Mach number distribution along the length of the pipe Is 
obtained, the other flow properties can be calculated by a simple rela- 
tion of the properties at two stations (1 and 2): 


Pi Ms / T 


k! 

y y - 

-»!) 

rti, 


1 


Zg = /T 


P P T 

Pi P1T2 ’ 


ii ■ i„ 


^“0 




- £2 

P.- ” Pi 

tl ^ 


y - 1 

2 ‘"2 -ry-l 


1 + M^ 


Z 2 -|> 

-1 + M? - 


The Mach number distribution In general depends on the nature of the 
temperature distribution and the heat flux along the pipe. For some 
cases, one can combine the friction coefficient and heat transfer coef- 
ficient via the Reynolds analogy, and the Integration of equation (56) 
might be possible In closed form [8], 
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D. Sudden Enlargement In a Pipe 


The only loss coefficient which can be calculated by simple 
analytical methods Is the case where there Is a sudden enlargement In a 
pipe. The problem was first treated by Borda and Carnot and Is there- 
fore known as the "Borda-Carnot loss." 

The flow Is again adiabatic, and It Is assumed that the pressure 
across the face of the enlargement Is equal to the pressure In the smaller 
pipe, just before the enlargement.* According to Figure lb, the three 
basic equations can be written as follows: 

Momentum 


m(v2 - Vi) = AgCP^ - Pg). 


(63) 


Continuity 


PxViAi = p2VaA2- 


(64) 


Energy 


ll + y = ll + — 21- £2 = 7 + 1 
2 7 ■ 1 Pi 2 7 ■ 1 P2 2(7-1) 


(65) 


where c* Is the velocity of sound where the Mach number Is unity. For 
adiabatic flow, the stagnation temperature Is constant, and we can write 
for V = c*, the right-hand term of equation (65) as 


Pi , 7 - , 1 . 
Pi 7 



y + 1 
2(7 - 1) 


Jz = y + ^ c*2 
2 27 



( 66 ) 


and 


Ps _ 7 + 1 *2 

P 2 27 



( 67 ) 


This assumption restricts the subsequent equations and M^ S 1. 
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From equation (63) and (64), we obtain 


~ _ Pa Pi 

PaVa^a Pa^a PiVi0* 


Hence, with equation (52) and (67) Introduced Into ( 68 ) 


^1 ~ "^a ^ ^ fy + 1 _ 7-1 1 _ fy + 1 _ 7-1 1 

c* V 2 \ 27 ” 27 ^*aj “ i 2 ^i \ 27 "27 ^*aj ’ 


where 0 is equal to the ratio Ai/Aa- In terms of the critical Mach 
number M*, we obtain 


( 68 ) 


(69) 



± fy + 1 

M* 1 27 

a ^ 


y - 1 

27 





(70) 


and 




0 . 


(71) 


This equation Is solved for M*> the Mach number In the larger pipe, and 
the positive sign of the square root Is taken. The density and pressure 
ratios can then be calculated as 


density - 



and, with the gas law (13), 


(72) 


pressure 




Z 

Z 

Z 

7 


- 1 
+ 1 
- 1 
+ 1 



(73) 
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The relationship between Mach number M and M*, is 


M^ = 




(7 + 1) - (7 - 1) M*‘ 


( 74 ) 


The density and pressure can be written as follows: 


density - 



0^1 / (7 - 1) Mg +~2 ~ 

^ Ma N/ (7 - 1) + 2 * 


(75) 


pressure 



m; 


Pi 

P2 


(76) 


E. A Sudden Contraction in a Pipe 

Here, the main features are an acceleration zone leading to the 
development of a vena contracta, followed by a deceleration zone similar 
to that analyzed by Borda and Carnot. 

The total pressure loss can be calculated approximately by the 
following idealization. The smaller tube with the area A2 and a given 
mass flow ifi is placed into a stream of velocity Vi < V2. Then the stagna- 
tion point of the free streamline is located on the outside of the tube. 
The gas contained inside the free stream tube of area Ai flows with an 
energy loss through the tube, since the flow separates at the lip of the 
pipe and contracts at the entrance (see Figure 10, lower part of the 
schematic). The surface of the free stream tube can now be solidified, 
thus forming the upstream tube with the area A^; (Ax > A2) . 

The momentum equation in the x-direction yields 


m V 2 + (P - Pi)2 A2 = m Vi + (P - Pi)i Ai + S + J (P - P^) dS^, 

S 

(77) 
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If we assume that the wall of the narrow pipe Is thin S/dg 0, 
then the suction force S can be taken as zero. The pressure difference 
(P - Pi) Integrated over the free stream tube surface Is zero by defi- 
nition of the free streamline. However, In practical cases, where the 
Idealization Is rather far off, a pressure difference (P - Pi), dependent 
on length b, builds up In the corner of the two pipes BCDE. With decreas- 
ing length b, the difference In the pressures P and Pi becomes more pro- 
nounced and the Integral must be taken Into account. One approximation 
could be obtained for the case of b > 0 and Ai » 

Abbreviating equation (77), we get 


T = 


S + 



(P - Pi) dS^. 


The mass flow through the tube Is constant and 


M, 


m* P 2 V 2 A 2 “ 7^2 (c /c )2 A 2 « 
®t 


(78) 


With the continuity equation. 


P^ 2 (c^/c )2 A2 “ PiMi(c^/c)i Ai, 


(79) 


and the relation. 


Pt2 _ 


. 12 

Pti P- ’ 


Pa/Pt2 -X 


(80) 


one can obtain, after rearranging equation (77) and solving for 


p 


1 1 +IQ- M^i(c^/c2)/(c^/ci) - T/(P2A2) 


the expression for the total pressure loss across the contraction 
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2 a 


1 + yK 


2 ^2. 

2 ^1, 


T 

PaAa 


The term T/PgAa is definitely a function of m| and Mf. For the case of 
Ml = 0, a very simple expression is obtained for calculating the Mach 
number Ma. We set 


^2 y Pti Pti 


p = 7 1 - 


Pa Aa 7 m: 


Then the Mach number in the pipe entrance becomes 


M? = 


(lea - 1' 


2 (7 - 1 - 2p2a) - (7 - 1 - (7 - 1 - 2p^a)2 ’ 


where the Ma is chosen, which is positive and smaller than unity. In 
engineering literature where energy loss coefficients are considered, the 
definition is 

, . 2^ ^ Pi - Pg . v! - vi = ^tl ~ ^t2 

^ ^ V| qa 

and in terms of the total pressure ratio between the stations 1 and 2 in 
Figure 1C, 


P - P P 

f* — ^ ^ ^ /p /p \ •2a c\ - 

I P 2 m| ^ 2 ^ 


where 


- (1 + M?) ; qa = ^ Pa 
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When we develop ^ series for small Mg and let Ai Increase 

Indefinitely, so that " 0, we obtain a simple expression for the 
energy loss coefficient 


5 ■ (i ■ - rMlp . K , fk - “'“=0 ■ 

Our next task Is to obtain an expression for the Integral term T. For a 
sharp-lipped Inlet of the pipe, the suction force S vanishes, and thus 

T (P - Pi) dS^. 

S 


In general, the pressure difference along the distance C - D In Figure 1C 
Is a function of and Mg. In our case, however, » 0, and therefore 
the term 


2T 

7 P2 As 



Cp(r/R2>(r/R2) d(r/R2) 


is only a function of the Mach number M 2 in the pipe. As a next step, 
we want to obtain an approximate expression for the pressure coefficient 
along the wall C - D. As a first approximation, which is valid for small 
Mach numbers M 2 , we treat the flow as incompressible, pi “ Pa “ p* The 
length b is zero. The Bernoulli equation along the streamline C - D 
yields 



2 P2 2 p 




and the pressure coefficient becomes 



For Ai ^00 , the velocity vanishes, and the pressure coefficient reduces 
to 
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The mass flow is constant between stations 1 and C, where Ag ■ KqA 2 and 
Kq is the contraction coefficient 


ih - p Ko As - p Va As 


and therefore 



The limit of the pressure coefficient between point C, which has shifted 
with Ai -» », and D is at point C V « “ 0, and the pressure coefficient 

vanishes. At point D, V « Vc = Vg/Kg and the pressure coefficient becomes 
Cp » -1 /Kq. Between the points C and D, the velocity behaves, at large 
distances from D, as a point source in the presence of a wall, and the 
velocity is proportional to l/Cr/R^)^. In the vicinity of D, r ■ R 2 , 
the velocity vanishes faster than 1/ (r/Rg)^, since the derivative of V 
at Ra is infinite; V, however, is finite and equal to Vq, The velocity 
derivative dV/dS, as a matter of fact, has theoretically a logarithmic 
discontinuity. If we still apply an exponential law to 


v(r/Ra) 


Vc 


1 


9 


the exponent n should grow indefinitely for r R^. 

If we accept the value of ^ « .5 for a sharp-edged orifice and 
for the length b = 0, as all engineering handbooks propose, an average 
value for 


n 


2 + K® 
o 


K 


s 

0 


will be obtained. For a vena contracta in a confined space, the pressure 
in the separation bubble is lower than in the atmosphere, tending to 
increase the contraction coefficient, Kq. According to figure 2, a value 
of Kq .8 is more appropriate than Kq .62, which is the value for a 
jet issuing from a flat-plate orifice into the atmosphere. With the value 
of Kq « .8, the exponent becomes n = 4.125, which seems to be reasonable. 
We choose n = 4.0, and the expression for the integral term becomes 


= Z i /[2 

P2A2 2 2 3k2 * 

o 
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III. THE ANALYTIC REPRESENTATION OF THE DISCHARGE COEFFICIENT 
AND OTHER FLOW AND ATMOSPHERIC PROPERTIES 


An analytical Investigation of the rather complex flow problems 
requires that the various coefficients of discharge, pipe friction, and 
the atmospheric properties, such as pressure and density, as a function 
of altitude, be represented In closed form, or given as tables. For 
efficiency, a closed analytical form of representation Is favored. 


A. Discharge Coefficient 

The discharge coefficient deviated remarkably from unity pri- 
marily because of the contraction of the streamlines. The streamlines 
curve or converge as they approach the orifice opening, bend around the 
edge of the hole and continue to curve and converge beyond the orifice. 

At some distance from the plane of the orifice, the het has a minimum 
section at which the streamlines are parallel. Because of the friction, 
the actual jet velocity in the orifice is less than the Ideal jet veloc- 
ity based on Isentroplc flow. A combination of these two factors Is the 
discharge coefficient. The coefficient increases substantially as the 
result of compressibility as shown in Figure 2. If orifices are built 
Into a pipe, the discharge coefficient depends also on the diameter ratio 
of the orifice and upstream pipe [5]. With orif ice-to-pipe diameter 
approaching one, the constriction of the flow disappears, and the dis- 
charge coefficient will also attain the value of one (Figure 3). 

Another important parameter is the Reynolds number formed by 
the orifice diameter and approach velocity. For small Reynolds numbers 
below 200, the discharge coefficient K decreases markedly as can be seen 
in Figures 4 and 5. Leak flows fall definitively into this region. 

Some of the curves In Figure 2 represent the free jet; l.e., 
the jet discharging into still air. (These curves are obtained from 
References 1 and 3.) A dependence on the aspect ratio for square open- 
ings is noticeable. For comparison, the discharge coefficient of a 
circular sharp edge orifice is also drawn. For the pressure ratio, 

Pex/F» approaching one, the coefficient has approximately the value of 
the theoretical calculations of the vena contracta of about « .6 - 
.61 by Treffts and .58 by Garabedian - [5]. The coefficients for two- 
and three-dimensional flow are probably the same. However, since the 
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combined effect of friction and contraction Is plotted, the difference 
can be attributed mostly to the first effect. The compressible jets 
have a higher curvature at the orifice [11], thus Increasing the minimum 
cross section of the jet far downstream of the orifice. 

For the subprogram, the discharge coefficient Kq was approxi- 
mated as 


K_ 


1 . 


.155 exp i -3,5(Pgjj/P^) 


.5186 


exp(Pg^/Pp + exp(-Pg^/pp ’ 

(84) 


which gives a fair value over the range 0 S ^ex^Pf ~ 

Figure 3 shows the change of the discharge coefficient obtained 
from Reference 5. At d/D = 1, the discharge coefficient Is Kq ■ 1, 
decreasing monotonlcally with d/D as predicted by theory to the value of 
Kq ■ .6, the discharge coefficient of an orifice cut into an Infinite 
wall. The theoretical value is obtained at approximately d/D ■ .23; 
then gravity effects become probably more pronounced for liquid jets. 

A fair approximation is given by the expression 


K = .6 + .4 exp 



(85) 


In Figure 4 the discharge coef flclentsfor plate orifices are 
drawn [10] , showing the influence of the Reynolds number for various 
orif ice-to-pipe diameter ratios. The Reynolds number is formed with 
the diameter and the approach velocity. 

The discharge coefficients for close-clearance orifices are 
drawn in Figure 5 for various Reynolds numbers obtained from Reference 6. 
The Reynolds number is here defined as 


R 


e 



A|i 




( 86 ) 
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where the hydraulic radius is given by 




area of flowing fluid 2 
wetted perimeter 


For a radial clearance, C, over flow length, L, approaching zero, the 
discharge coefficient for all Reynolds numbers decrease. For c/L ^ 1, 
the coefficient is only a function of the Reynolds number. 

A general approximate expression for K/Kg as a function of C/L 
and Rg (Figure 5) is given by 


K /<, -4.788 ^/c7L^ , ^ 1 

— “ (^1 - e ^ (a + n In Rg). 


(87) 


The constants differ for different Reynolds numbers and are given by: 


1 < Rg ^ 100 a = .02 

100 < Re S 500 a = .52 

Rg > 500 a = 1.0 


n = .1867 
n = .0782 
n = 0. 


(The hump in the curves at low C/L and high Rg numbers is not considered 
In these expressions.) 

The coefficient can be set equal to .62 without considering 
any effect of the outside flow. 

When the jet Issues into the outside stream, flowing normally 
to the jet axis, afurtherdrop of the discharge coefficient is expected. 
In Figure 6 the discharge coefficient ratio, K/K^, is drawn for an 
orifice [1] of aspect ratio 1 and for three different free stream Mach 
numbers over the mass flow rate ratio (m; Kq/hIooK). A definite influence 
of the free stream Mach number M^ can be observed. In Reference 1 the 
test was conducted with orifices of different aspect ratios and shapes 
cut into a wall placed tangentially to the outer flow field. At the 
position of the orifice, the boundary layer was relatively thin. With 
thicker boundary layer flow regions in the vicinity of the orifice, the 
discharge coefficient K should be somewhat improved. Further tests are 
necessary in this area. 
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For the venting problem, a mean representative curve was derived. 
Its expression is 


_K 

K. 


“i^o 

Cl - X'>E + X*' ’ ^ ” 00 ^ P. 


X‘ 



( 88 ) 


where n * 4.5714 and E = .0101 representing the curve for M = 1 closely. 

To calculate the proper discharge coefficient, it is necessary 
to determine the absolute viscosity of the discharging gas. In general, 
the viscosity, n, is a function of the pressure and temperature. How- 
ever, as the pressure and thus the density of the gas decreases, the 
absolute viscosity approaches the low density limit and the pressure 
dependence becomes less. For most gases, this low density limit is 
reached at a pressure of about one atmosphere. Since the pressure in 
the compartment is about one atmosphere at the beginning of venting and 
decreases further with time, we consider therefore only the temperature 
dependence. 

The dependence of the viscosity on the temperature, obtained 
from Reference 4, is plotted in Figure 7 as a mineral plot of all gases. 
It gives the relationship between the reduced viscosity, |ij^ = |i/|ic» 
reduced temperature, Tj^ = T/Tg with the reduced pressure P,. = P/Pg as 
parameter; Hg, Tg, and Pc are the critical viscosity, temperature, and 
pressure, respectively, of the particular gas under consideration. 

Experimental values of Hg are seldom available. However, (jg 
may be estimated in the following way; 


IJi^ = 7.7 P§/^ Tg^/® X lO'"^ [Kg/m sec], (89) 


where Pg is the critical pressure in atm (Kg/ cm^) and Tg is the critical 
temperature in °K. 

The dependence of on T^ is in general 

- (a V". (90) 
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where the constants are for the regions 


2.0 < S 

3.0 < S 

5.0 < T_ ^ 

K 


1.5 

2.0 

3.0 

5.0 
10.0 


a - .442 
a “ .462 
a = .500 
a - .564 
a - .631 


n - .9882 
n - .8204 
n - .6480 
n ■ .5185 
n * .5099. 


B. The Atmospheric Properties 


For the calculation of the atmospheric properties of air, the 
formula and constants of Reference 9 were used. However, the rather com- 
plex expression for the acceleration due to gravity was further approxi- 
mated by 


g 


®o [a + z]^ 


(91) 


where a = 6,378,178 meters is the radius of the earth for the geographic 
latitude 0 = 45® and z is the height above the earth's surface. The 
geopotential altitude becomes then 



(92) 


Integration yields 


H 


a 


z 

a + z * 


(93) 


A comparison of the values obtained by this rough formula with the values 
of Reference 9 showed an error of 3 meters for an altitude of z = 70,000 
meters. Other properties are obtained by the following formulas: 
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Pressure 




g H 
00 


mb 

+ L' (H - H, 






for l 4 ^ 0 


(94) 


or 


P ■ exp 


Sc “o"-' 


R*T 


mb J 


for “ 0. 


(95) 


Density 


M P 
o 



(96) 


Temperature 


^ ■ ’■mb + ’■> - «b>- 

Speed of Sound 



(97) 


(98) 


In these equations the subscript b designates the respective values at 
the base of the sublayer. The constants for all layers are 


= 9.80665 m/sec^ 

Mq = 28.9644 Kg/K mol 

R* = 8.31432 Joules (°K)‘^/mol. 
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The constants for the different sublayers are 


H = 0 + 11 Km 

- -6.5; - 288.15*K; - On; P,, - 1.013250 x lO^Newton/m^ 

H = 11 + 20 Km 

L' = 0; T . - 216.65“K; 
m mo 

H = 20 -5- 32 Km 

L' = 1; T . - 216.65®K; 

H = 32 -s- 47 Km 

= 2.8; = 228.65“K; = 32000 m; * 8.68014 mb 

H = 47 -5- 52 Km 

^m “ ^mb “ 270.65®K; Hj^ = 47000 m; * 1.10905 mb 
H ® 52 + 61 Km 

= -2; = 270. 65 “K; = 52000 m; P^^ = 5.90005 mb, etc. 


- 11000 m; P^^ » 2.26320 x 10® mb 


- 20000 m; P^ » 5.47487 x 10 mb 


C. The Friction Factor 


In the equations for pipe flow, the friction factor, f, appears; 
accurate numerical values for this factor will be needed for the solution 
of engineering problems which require information on energy loss. Notice 
that f is an experimental coefficient usually determined by energy losses, 
the length and diameter of the pipe, and the velocity. 

Nicuradse [10] conducted systematic tests, using measurable 
roughness produced by uniform sand grains of diameter e. Through the 
use of such artificial roughness, he was able to show that the friction 
factor depended upon Reynolds number and the relative roughness e/d (Fig- 
ure 8). In this diagram the friction factor, f, is already the total 
value for a pipe of length, L. The Reynolds number therefore uses the 
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total length, L. At high Reynolds number and high relative roughness, 
the friction factor becomes wholly a function of the relative roughness, 
e/d. Colebrooh and White [10] were able to correlate the results of 
many laboratory tests on commercial pipes over a wide range of Reynolds 
number and roughness by the following expression: 


^-nos.o|-l.U 


21ogio 


'l + 1. 

Rg(e/d) •JTT -I 


(99) 


However, care must be taken In defining relative roughness. For Instance, 
If a pipe has a definite "regular" roughness, like a wavy wall and other- 
wise a smooth surface, then the curve would show a pronounced tendency to 
vary with the Reynolds number and a friction factor for a hydraulically 
smooth pipe has to be taken. 

For a fully developed velocity profile, l.e., at distances 
from the pipe Inlet greater than 50 diameters, no significant effect of 
Mach number was observed. The pipe Reynolds number is given as 


Rg = vp D g/p. 


For short pipes near the duct inlet, the one -dimensional analysis is 
misleading, since the velocity profile is not yet fully developed and 
changes in momentum flux of appreciable magnitude may occur. At the 
entrance, the flow is rather similar to the flat-plate flow, and the 
Reynolds number based on the flow length from the inlet is more important. 
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FIG. 2. DISCHARGE COEFFICIENT 
AS A FUNCTION OF COMPRESSIBILITY 


CoNtrcetlM CMfficitiit, k 



Orlfie* Dit«/PiM DiiM, 4/0 

FIG. 3. DISCHARGE COEFFICIENT AS A FUNCTION OF 
ORIFICE-TO-PIPE DIAMETER 
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FIG. 4. DISCHARGE COEFFICIENTS FOR THE PLATE ORIFICE 


radial clearance 


flow length L 
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( rodial clioronci )/( flow longth ) 

FIG. 5. LEAK DISCHARGE COEFFICIENT AS A FUNCTION OF 
REYNOLDS NUMBER AND RADI AL CLEARANCE TO FLOW LENGTH 
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Roughneii of pipe ii meoiured by <• end boi 
typical values as folloKS (after Moody): 


Pipe 

c.ft 

Drowp tubing 

.000005 

Commercial steel 

.00015 

Asphalted cast Iron 

.0004 

Galvanized Iron 

.0005 

Cast Iron 

.00085 

Concrete 

.00I-.01 

Riveted steel 

.003-.03 


FIG. 8. FRICTION COEFFICIENT VERSUS 
PIPE REYNOLDS NUMBER FOR 
INCOMPRESSIBLE, FULLY DEVELOPED FLOW 
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TABLE 1 

CRITICAL PROPERTIES 


Substance 

Molecular 

Weight 

M 



Critical Constants^*®*’* 


T. 

(”K) 

(atm) 

(cm’ g-molc"‘) 

P4 

(g cm~‘ scc'O 

X 10* 

A*. 

(cal see”* cm”* * K-») 
X 10* 

Lijiht tU’menu! 
H, 

2.010 

33.3 

12.80 

65.0 

34.7 


lie 

4.003 

5.26 

2.26 

57.8 

25.4 

— 

Noble gases: 
Nc 

20.183 

44.5 

26.9 

41.7 

156. 

79.2 

Ar 

39.944 

151. 

48.0 

152 

264. 

71.0 

Kr 

83.80 

209.4 

54.3 

92.2 

396. 

49.4 

Xc 

131.3 

28^.8 

58.0 

118.8 

490. 

40.2 

Simple polyatomic 
substances: 

Air 

28.97 

132. 

36.4 

86.6 

193. 

90.8 

N, 

28.02 

126.2 

33.5 

90.1 

ISO. 

86.8 

o. 

32.00 

154.4 

49.7 

74.4 

250. 

105.3 

o. 

48.00 

268. 

67. 

89.4 

— 



CO 

28.01 

133. 

34.5 

93.1 

190. 

86.5 

CO, 

44.01 

3at.2 

72.9- 

94.0 

343. 

122. 

NO 

30.01 

ISO. 

64. 

57. 

258. 

118.2 

N,0 

44.02 

309.7 

71.7 

96.3 

332. 

13J. 

SO., 

64.07 

430.7 

77.8 

122. 

411. 

98.6 

F. 

38.00 

— 

— 

— 

— 



Cl, 

70.91 

417. 

76.1 

124. 

420. 

97.0 

Br, 

159.83 

584. 

102. 

144. 

— 

— 

I. 

253.82 

800. 

— 

— 

— 

— 

Hydrocarbons: 

C\h 

16.04 

190.7 

45.8 

99.3 

159. 

158.0 

Qii, 

26.04 

309.5 

61.6 

113. 

237. 

— 

QH 4 

28.05 

282.4 

50.0 

124. 

215. 

— 

QH. 

30.07 

305.4 

48.2 

148. 

210. 

203.0 

QH, 

42.08 

365.0 

45.5 

181. 

233. 

— 

C,K, 

44.09 

370.0 

42.0 

200. 

228. 

— 

rt-C|Hn 

58.12 

425.2 

37.5 

255. 

239. 

— 


58.12 

i 408.1 

36.0 

263. 

239. 

— • 

//•CjHji 

72.15 

469.8 

33.3 

311. 

238. 

— 

/;-QHu 

86.17 

507.9 

29.9 

368. 

248. 

— 


100.20 

540.2 

27.0 

426. 

254. 

— 

/j“C,Hji 

114.22 

569.4 

24.6 

485. 

259. 

— 

rt-CjH,® 

128.25 

595.0 

22.5 

543. 

265. 


Cyclohexane 

84.16 

553. 

40.0 

308. 

284. 


C.H. 

78.11 

562.6 

48.6 

260. 

312. 

— 

compounds: 

CH» 

16.04 

190.7 

45.8 

99.3 

159. 

*58.0 

CH.Cl 

50.49 

416.3 

65.9 

143. 

338. 

-- 

CH. Cl, 

84.94 

510. 

60. 

— 

— 

— 

CHCl, 

119.39 

536.6 

54. 

240. 

410. 

— 

CO. 

153.84 

556.4 

45.0 

276. 

413. 

-- 

QN, 

52.04 

400. 

59. 

— 

— 

— 

COS 

60.08 

378. 

61. 

— 

— 

— 

CS, 

76.14 

552. 

78. 

170. 

404. 

— 
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APPENDIX 

PREPARATION OF THE INPUT DATA FOR THE DIFFERENT VENTING PROGRAMS 


A. Venting of Compartments In Series 

The simple discharge of a number of compartments In series Is 
calculated with this program. The single compartments are connected by 
only one orifice. The first compartment, which vents into the atmosphere, 
can have more than one discharging orifice. For all compartments the 
possibility of leak flow is provided, with the restriction, however that 
only one leak per compartment can be handled. 

All input data are read in by the INPUT subroutine. Floating- 
point numbers have the Format E15.8 and fixed-point numbers have I3. The 
data in consecutive order have the following meaning and dimensions. 

NH, NP, MPRNT, NPO 

NH = Fixed-point number, designating the number of trajectory 
data cards. 

NP = Fixed-point number, designating the number of cards for 

the pressure data outside of the leaks of the compartments. 

MPRNT = Fixed-point number for printing out check data. MPRNT 
greater than zero for print out. 

NPO ■ Number of cards for the pressure data outside the 

orifices of the first compartment. The first compart- 
ment vents into the atmosphere. 


OM, TC, PC, GC, G 

CM = Molecular weight of the vented gas. 

TC = Critical temperature of the vented gas in [®Kj. 

PC = Critical pressure of the vented gas in atmospheres 

[Kg/cm^] . 

GC “ Gas constant of the vented gas [m/®K]. 

G = Specific heat ratio of vented gas. 
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TT(I), HT(I), FMT(I) NH times, maximum NH - 40 

TT(I) - Trajectory time in seconds. 

HT(I) - Trajectory altitude in meters m at time TT(I). 

FMT(I) ■ Trajectory Mach number at time TT(I). 

FM(J), CPF(J,1), CPF(J,2), ... NP times, maximum NP “ 40 

FMF(J) = Mach number for which the pressure coefficients 
CPF(J,I) are selected. It is not necessary that 
this Mach number correspond to the FMT(I) mentioned 
above. Both, however, must be in the same range. 

It is recommended that more points be selected in 
the transonic range. 

CPF(J,I) = Pressure coefficients of the leaks, outside the com- 
partment. If more compartments are to be vented, 
other cards must be added and the dimensions in the 
COMMON statements must be changed accordingly. 


FMO(I), CP0(I,1), CP0(I,2), ... NPO times, maximum NPO = 40 

FMO(I) = Mach number for which the orifice pressure coefficients 
C0P(I,J) are selected. 

CP0(I,J) = Orifice pressure coefficients of the first compartment 
outside the compartment. 


NO, NT, N, NOR 

NO = Fixed-point number, designating the first time step 

minus 1 second. 

NT = Fixed-point number, designating the end time step 

minus 1 second. 

N = Number of compartments in series also total number of 

leaks; one leak per compartment. 

NOR = Fixed-point number; number of orifices of the first 

compartment venting into the atmosphere. 
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NOR times 


AO (I), G1(I), OMl(I) 

AO (I) ~ Area of outside orifices of the first compartment 

[m®] . 

Gl(I) Mass flow through orifice number I of the first 

compartment [Kg sec/m] . 

OMl (I) = Mach number of orifice flow. 

For the very first time step, insert a zero for Gl(I) and (Ml(I). After 
NT time steps, the OUTPUT subroutine will punch all cards beginning with 
the fixed-point numbers NO, NT, N, NOR and later. 


A(I), VO(I), R1(I), PI (I), AMI (I) 

AL(I), HR (I), CCF(I) N times 

A (I) = Orifice area of the orifices connecting the compart- 

ments [m^] . 

VO(I) = Volume of the i-th compartment [m®]. 

Rl(I) = Initial density of the gas in the compartments 

[Kg sec^/m'^]. 

Pi (I) = Initial pressure of the gas in the compartments in 

[kg/m^] . 

AMI (I) ■ Initial mass flow out of the compartment [Kg sec/mj. 

AL(I) = Leak area of the compartments (representative mean 

value) [m^]. 

HR(I) ■ Hydraulic radius or mean radius of leaks defined in 
III. A [m®]. 

CCF(I) “ Radial clearance over flow length of the leaks 
according to Figure 5. 

The program can be restarted at any time step NT. The punched data 
replace the old corresponding data. The program punches for NT a fixed- 
point number. It can be changed to any other value NT ^ NO. 
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B. The Input Data for the Venting Program with Pipe Flow 


In this program, the pipe replaces the first compartment. Some 
of the Input data have been changed. The format statements remain the 
same as In A. We repeat here also those symbols which have not changed. 


NH, NP, MPRNT 

Fixed-point numbers with the same meaning as In A. 

(M, TC, PC, GC, G 

Same meaning as In A. 

TT(I), HT(I), PMT(I), QT(I) NH times, maximum NH = 40 

Same as In A. For the program with pipe flow and heat addi- 
tion the term QT(I) Is added. 

QT(I) ■ Total heat flux per unit mass of gas [Kg m/sec]. 

FMF(I), CPF(I,1), CPF(I,2), . . . NP times, maximum NP = 40 

Same as In A. 

NO, NT, N 

Same as in A. 

A(l), HR(1), PL(1), Rl(l), Pl(l) 

A(l) = Cross section of the pipe [m^] . 

HR(1) = Hydraulic radius of pipe [m] . 

PL(1) = Pipe length [m] , 

Rl(l) = Initial density of the gas in the pipe [Kg sec^/m"*^]. 
Pl(l) = Initial pressure of the gas in the pipe [Kg/m^]. 
These data have been added to the original program. 
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A(I), VO(I), R1(I), PI (I), AMI (I) 
AL(I), HR(I), CCF(I) 


N-1 times 


Same as In A. 


W 

W additional loss coefficient of the pipe as bends and elbows, 
etc. If no additional losses, enter zero. 


C. Description of the Iteration Subroutine 

This subroutine determines the next Iterated pressure at the 
entrance of the pipe joining the upstream compartment. For subsonic 
exit Mach number the pressure at the exit has to match the external pres- 
sure. For this case the entrance pressure Is varied according to a cer- 
tain Iteration scheme until the condition of equal pressures at the exit 
is met. For supersonic pipe flow the exit is choked and the pipe exit 
pressure is higher than the external pressure. Further expansion of the 
flow takes place external to the pipe. 

In the subroutine the first trial value of the entrance pressure, 
which is determined in the main program by an extrapolation of entrance 
pressures occurlng at an earlier time step. This pressure causes a 
certain exit pressure and exit Mach number, calculated by the subprogram 
PIPE. This program determines also an entrance Mach number which corre- 
sponds to a choked exit. The pressure difference across the exit 
* Pgjjit " ^exlt 00 Mach number difference at the entrance 

AM= Mgj^^j. - Mentr^^ex “ listed according to whether Mgj^ is smaller 

than one (OK = AP) or equal to one (OG = AP) . To obtain the second value 
of Pgntr* first one is perturbated by a fraction of the pressure 

difference " ^entr^ ’ schematic Al. AP and M is again listed. The 

next improved value is obtained by extrapolating the first and second 
approximation to make at the entrance of the pipe equal to zero. If 
zil has reached the tolerance, a logical decision is made. For AM ^ Tol 
and AP < 0 the pipe exit is not choked and the condition AP = 0 must be 
sought. This condition is met by extrapolating PK (for which Mgj^ < 1) 
to OK = AP = 0. If, however, M ^ Tol and P ^ 0, the condition of choked 
exit flow is met within the tolerance. 
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pH) p(*)\ 


SCHEMATIC A-1. SOLUTION FOR 
THE TWO DIFFERENT PIPE FLOW REGIMES f'Af „jt ^l) 





PROGRAM VENT 

c venting of compartments in series 

C «***»**««*«««*««*««*«'»«««««*«*«*««*««««« 

C NOR SNUMBER OF ORIFICE OF THE FIRST COMPARTMENT 

C NPO » NUMBER OF CARDS OF PRESSURE DATA OUTSIDE OF ORIFICE 

C NH = NUMBER OF CARDS FOR TRAJECTORY DATA 

C N -NUMBER OF COMPARTMENTS » ALSO NUMBER OF LEAKS 

C ONE LEAK PER COMPARTMENT 

C PIP^E IS CONSIDERED AS COMPARTMENT NO* I 

C HT - ALTITUDE AT TIME TT IN METERS 

C FMT - MACH NUMBER AT TIME TT 

C A - ORIFICE AREA BETWEEN N COMP. IN METERS 

C VO = VOLUME OF COMPARTMENTS IN METERS3 

C CPF- PRESSURE COEFFICIENT AT MACH N. 

C A « AREA OF NV ORIFICES OF COMP. 1 IN METERS 
CP- PRESSURE IN COMPARTMENT IN (KG/M2) 

C CM1= COMPARTMENT MACH NUMBER 

C ALl- LEAK MASS FLOW 

C R = DENSITY OF N COMP. IN KG.'SEC./MA 

C AL -LEAK AREA OF COMP.N IN METERS 
C HR -hydraulic RADIUS OR MEAN RADIUS IN METER 

C TC = CRITICAL TEMPERATURE OF GAS 

C GC = GAS CONSTANT OF GAS (M/DEG. K) 

C PC = CRITICAL PRESSURE OF GAS IN ATM(KG/CM2) 

C QM » MOL WEIGHT OF PARTICULAR GAS 
C CCF- RADIAL CLEARANCE OVER FLOW LENGTH 

C WI -ADDITIONAL LOSS COEFFICIENT FOR PIPE* IF NONE ENTER 0. 

C IN PROPER FORMAT 

C NUMBER 1 COMP. IS THE COMP. THAT VENTS INTO THE ATMOSPHERE 

C AMI = MASS FLOW OUT OF THE COMPARTMENT 

C G1 = MASS FLOW OUT OF ORIFICES OF THE FIRST COMPARTMENT 

C PG = PRESSURE OUTSIDE OF THE ORIFICES OF THE FIRST COMPARTMENT 

COMMON MPRNT*N0R»NP0 

COMMON NH»N»NT » NP » GC »0M * TC ♦ PC *G» TT ( AO ) »HT ( AO ) » FMT (AO) *FMF(AO) *W»NO 
COMMON CPF(A0»5) »A(5) .HR(5) »PL(1) »V0(5) »AL(S) »CCF(5) »PAl(t>) »P2( i) 
COMMON P1(5)»R1(5) »AM1(5) .T1.CMK5) *ALK5) »G1 ( 5 ) *0M1 ( 5 ) 

COMMON P2(5).R2(5).AM2(5)»T2»CM2(S)*AL2(S)*G2(5)*OM2(5) 

COMMON PG(5) »CPO(AO»5) »FMO(AO) »AO(5) *PA » TA *RA »CA* ALT ♦ FM » WH( AO * 5 ) 
COMMON TP (AO) »OP(AO) »WM(A0*5) 

DO 5 I-1.5 
AM1(I)=0.0 
5 AK2( I )=1.0E~07 
CALL INPUT 

D=(2./(G+1. ) )**(G/(G-1. ) ) 

DO 13 I-NO.NT 
11=1-1 
T 1-1.1 
T2-T1-1.0 
CALL TRAJ 
IF( I-1)A»11»A 
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4 DO 15 L=l»40 

IF(MPRNT-4)24»24»25 
25 WRITE(51»555) 

WRITE(61#90)L»MPRNT 
24 IF{L-1)6*6,7 

6 00 8 K>1*N 

8 AMI (K»»AM2(K) 

7 DO 20 K»1*N 
IFU-1>1»1»2 

1 CONTINUE 

CALL PARTS! I .K*L»D) 

GO TO 23 

2 CONTINUE 

CALL PARTA! I *K.L.D) 

23 WH(L»K)»P.i(<-l) 

20 WM(L»K)=AM1 (K) 

IF( L-1 ) 10»10»9 

9 DO 3 K=1»N 
IF(MPRNT) 18»18» 19 

19 WRITE(61*91 )K»WM(L*K) .WM(L-1»K) 

18 IF(ABS(WM(L*K)-WM(L-1»K) )-l*0E-07) 3 *3»10 

3 CONTINUE 
GO TO 12 

10 IF(MPRNT)21»21*22 

22 ,WRITE(61»90)L»L#AM1! 1) »AM1(2) .AMKO) *AM1(4) 

21 IF(L-3) 15»16.16 

16 DO 17 K=3.N 

17 PKK.-1 ) *WH( L.K)-(WM(L.ia-WM(L-l»K) ) * ( WH ( L » K) -WH ( L“1 »lO )/(WM(L»K ) 
12iO*WM(L-l.K)+WM(L-2»K) ) 

15 CONTINUE 
GO TO 12 

11 DO 14 IP=1 *N 
CMK IP)=+0.0 

14 ALl(IP)=+0.0 

12 CALL OUTPUT ( I ) 

WRITE«61»555) 

13 CONTINUE 
STOP 

555 FORMAT! IHO) 

90 FORMAT(2I4.8E15.8) 

91 FORMAT! I4*5E15. 8) 

END 
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c 


SUBROUTI NE PARTAd *K«L»D) 
VENTING OF COMPARTMENTS IN SERIES 
DIMENSION T(40)*O(A0) 


COMMON MPRNT»NOR*NPO 

COMMON NH»N»NT*NP»GC.OM»TC»PC.G»TT(40) »HT ( 40 ) »FMT ( 40 ) »FMF ( 40 ) ♦ W*NO 
COMMON CPF(40»5) »A(5) »HR(5) .PLU) tVOIS) »AL(5) #CCF<5) ♦PAK5) iP^I D' 
COMMON Rl(5) »R1 (5> tAMKS) »T1.CM1(S) .ALl ( 5 f *61 ( 5 ) .OMl ( 5 ) 

COMMON P2(5) *R2(5) »AM2{5) *T2*CM2(S) *AL2(3) *G2 ( 9 ) *OM2 ( S ) 

COMMON PG( 5) »CPO(40*9> *FM0(40) *AO( 9) »PA*TA*RA*CA*ALT»FMtWH(40«5 ) 
Tm**P2(K)*.999 


K1»K+1 


DO 30 Mal»40 

A1=(P1(K-1)+P2(K-1) )/2«0 
RKK)=R2(K)*( (T(M)/P2(K) )**<1./G)) 

CALL APPT(K»A1»R1(K) »T(M) .A (K>»A4»D) 

IF(AL(ia )27»27.33 

27 AL1(K)=0.0 
GO TO 28 

33 21*PA/T(M) 

CALL LEAK(21»T(M) ♦RIU) #AL(K) *0»HR I K ) »CCF( K ) . ALl ( K ) ) 

28 Sl=R2(K)*VO(K) 

S2=(2.0*( T(M)/P2(K) )-1.0 )**( l./G) 

S3»S1*( 1.0-S2 )+AMl(Kl )-ALKK> 

CMX2=2.0*< (T(M> /Al)**( ( G-1 . ) /G ) -1 . 0 ) / ( G-1 . ) 

CMl (K)«SQRT (ABS(CMX2) I 
IF(CM1(K)-1.0)6»6»7 
7 CM1(K)=*1.0 

6 0(M)=A4-S3 

IFIM-l )54»54»35 
54 T(M+1)=TIM)«1.005 

GO TO 62 

35 T(M+1 )aT(M)-0(M)*(T(M)-T(M-l) ) / ( 0( M ) -0(M-1 ) ) 

62 PI (K)=2.0*T (M+1 )-P2(K» 

AM1(K»=S3 

IFIABS (0<M) )-1.0E-07)22»30»30 
30 CONTINUE 

WRITE(61«99) 

22 IF(MPRNTJ4»4.5 

5 WRITE (61 *80) I »K « L*M«D»T (M ) *P1 ( K ) *S3 »0(M) »R1 ( K ) 

4 RETURN 

80 FORMAT(4I4,6E15.8) 

99 F0RMAT(17H PARTA TOLERANCE) 

END 
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SUBROUTINE PARTB ( I >K >L >D ) 

C VENTING OF COMPARTMENTS IN SERIES 

DIMENSION T(A0) »0(40) #DO<40) *TO(AO) 

COMMON MPRNT*NOR*NPO 

COMMON NH*N»NT.NP.GC»OM»TC*PC»G*TT(40) »HT(40 J »FMT(40) »FMF ( 40 ) t W «NO 
COMMON CPF(40»5) ♦ A ( 5 ) *HR ( 5 ) »PL ( 1 ) »V0 ( 5 ) » AL ( 5 ) *CCF<5) » PAl < 5 ) ♦PH ( 1 J 
COMMON Pl<5) »R1<5) *AM1<5) *T1 *CM1 ( 5 ) »AU ( 5 ) »G1 < 5 ) *0M1 ( 9 » 

COMMON P2( 5) «R2 (5) »AM2( 9 ) »T2»CM2(5 ) >AL2( 5) *G2 ( 5 ) *0M2 ( 9 ) 

COMMON PG(9) •CPOI40f9) iFMOI40) »A0(9) i TA *RA f CA *ALT * FM * WH ( 40 * 9 ) 
COMMON TP(40) f0P(40) »WM(40*5) 

T( i)=P2(K)*.999 
K1*K+1 

DO 30 M=l«40 

Rl(ia=R2lK)*( (T(M)/P2(ia )**a./G)) 

SUMa.O 

DO 6 J=1*N0R 

CALL APPT(J»PG( J) *R1 (lO »T(M) »A0( J) »A4*D) 

A1=PG( J)/T(M) 

CALL COEF(Al.AK) 

AOMsRA*FM*CA 
Y = AK 

DO 3 Jl=l*30 
TOC Jl)aY 
IF(Al-D) 1»1»2 

1 A1»D 

2 RI»R1CK)*(A1»*( l./G) ) 

AI=SQRT(ABS(6*T(M)*R1CK) ) )' 

VI=AI*SQRT( ABS(2«*( 1.tA1**( (G-1. ) /G)) )) 

AIM=RI*VI*AK 

X=AIM/(AOM*Y) 

CALL VENTCC X*AK»Y1 ) 

DOC J1 )=Y1-Y 
IFC Jl-2»8,8»9 

8 Y=.95*AK 
GO TO 12 

9 Y=TOC J1)-D0CJ1)*CT0CJ1)-T0CJ1-1) ) / C DO C J1 ) -DO C Jl-1 ) ) 

12 IFCMPRNT-4) 10*10»13 

13 WRITEC61.81)J.J1*T0CJ1).Y»D0CJ1).Y1»A1»D*X 

10 IFCABSCDOC Jl) ) -1 .OE-0 8 ) 1 1 . 1 1 . 3 

3 CONTINUE 

11 G1 C J) =A4*Y/AK 
SUM = SUM+G1 C J) 

CMX2 = 2.0»C C TCM)/PGC J) )**C C G-1 . ) /G ) -1 . 0 ) / C G-1 . ) 

OMIC J)=SQRTCABSCCMX2) ) 

IFCMPRNT-3) 14*14.15 

15 WRITEC61.80) J*M»K»J1»A1»AK»A4»0M1C J) *G1C J) »TCM) 

14 IFCOM1CJ)-1.0)6*6»7 

7 OM1CJ)=1.0 

6 CONTINUE 

IFCALCK) J27.27.33 
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27 ALl(K)-0.0 
GO TO 28 

33 21»PA/T(M) 

CALL LEAK(Z1*T(M) »R1(K) »AL(iC) *0*HR ( K) *CCF( K) » ALl ( K ) ) 

28 Sl-R2(K)*VO(iO 

S2« (2.0<M T(M)/P2(K) ) -1 »0 ) ** 1 1 . /G ) 

S3-S1*( 1.0-S2)+AM1(KX )-ALl(K) 

0<M)*SUM-S3 
IF(M->I)34»34*33 
54 T(M+1)-T(M)*1.005 

GO TO 62 

35 T(M+1)»T(M)-0(M)*(T(M)-T(M-1) )/(0(M)-0JM-l) ) 

62 Pl(K)«2.0*T(M+l)-P2U4 

AM1(K)°S3 

IFtABS (0(M) )-l*0E-07)22»30»30 
30 CONTINUE 

WRITE(61»99) 

22 IF(MPRNT)4»4»5 

5 WRITE<61*80) I *K*L*M»T(M) »Pl(iO *S3*0(M) »Rl(lO *SUM 

4 RETURN 

80 FORMAT(4I4»7E15.8 ) 

99 FORMAT(17H PARTA TOLERANCE) 

81 FORMAT(2I4»8E15.8) 

END 
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SUBROUTINE TRAJ 
pA»AMBIENT I^R'eSSURE 
TA=TEMPERATURE 
RA=DENSITY 
CA»SPEEO OF SOUND 

PAl (N>«PRESSURE OUTSIDE PIPE AND LEAKS AT TIME T1 
DIMENSION WK40) 

COMMON MPRNTtNOR*NPO 

common NH*N*NT»NP.GC*OM*TCfPC»0»TT(<fO) »HT(40» »FMT(40) »FMF(40J »W»N0 
COMMON CPF (40. 5 > . A ( 5 ) .HR ( 5 ) » PL ( 1 ) . VO ( 5 ) ♦ AL ( 5 ) .CCF ( 5 ) » PAl ( 5 ) » PZ ( 1 ) 
COMMON P1(5).R1(5).AM1(5) .TI.CMK5) .ALK5) .61(5) .0MK5) 

COMMON P2( 5) .R2 (5) .AM2(5) .T2.CM2 (5> »AL2( 5) .02 I 5 ) .0M2 ( 5 ) 

COMMON PG( 5) .CPO(40.5) .FM0(40) .A0(5) »PA»TA»RA»CA»ALT»FM 
CALL INTER! 1 . NH . 1 . 1 . TT .HT .T1 . ALT ) 

CALL INTER! 1 . NH . 1 . 1 . TT . FMT .T1 .FM) 

DO 1 I=1.N 
DO 2 J=1.NP 
W1!J)=CPF!J.I ) 

CALL INTER! 1 .NP . 1 . 1 . FMF . W1 .FM.PAl ( I ) ) 

CONTINUE 

HG=ALT*6.378178E+06/ (6.378178E+06+ALT) 

CALL ATMOS!HG.PA.TA.RA.CA) 

DO 3 I=>1.N 

PAl! I )=PA+RA*PA1< I )*( (FM*CA)#*2)/2.0 
DO 4 1=1. NOR 
DO 5 J=1.NP0 
Wl! J)=CPO!J.I ) 

CALL INTER! 1 .NPO. 1 . 1 . FMO. Wl .FM.PG! I ) ) 

CONTINUE 
DO 6 1 = 1. -NOR 

PG! I )=PA+RA*PG( I )*( ( FM*CA)**2 ) /2.0 

RETURN 

END 



SUBROUTINE INPUT 

C VENT I NG OF C6K<i'A'RTMENTS IN SERIES 

C NP » NUMBER OF CARDS FOR PRESSURE DATA OUTSIDE THE LEAKS 

C NH » NUMBER OF CARDS FOR TRAJECTORY DATA 

C NPO » NUMBER OF CARDS FOR PRESSURE AT ORIFICE 

C FMO » MACH NUMBER FOR PRESSURE DATA OUTSIDE THE ORIFICE 

C CPO = PRESSURE COEFFICIENT OUTSIDE THE ORIFICE 
C AO s ORIFICE AREA 

C PG a PRESSURE OUTSIDETHE ORIFICE 

COMMON MPRNT»NOR»NPO 

COMMON NH»N*NT»NP»GC»0M»TC»PC»G»TT( A0> »HT(40) »FMT(40) »FMFI40) »W*NO 
COMMON CPF(40.5) »A(5) .HR(5) *PL(1) ♦VO(5> fALIS) »CCF(5) »PA1(5) »P2( 1) 
COMMON PI ( 5) »R1 (5 ) »AM1 (5 ) ^TltCMl (5) »AU1 (5) »G1 ( 5) *OMl(5 ) 

COMMON P2(5» »R2t5) ♦AM2(5) »T2fCM2{5) ♦ AL2 ( 5 ) »G2 ( 5 ) *OM2 ( 5 ) 

COMMON P6( 5) »CPO(40t5) »FM0C40 ) »AO( 5 ) »PA » TA »RA » CA f ALT • FMTWm^O » 5 ) 

READ(60»100)NH»NP.MPRNT*NPO 

READ(60.101 )OM.TC.PC»GC*G 

VC=7*7*SQRT (OM)*(TC**«-l./6. ) )*(PC**(2./3. ) )*1.0E-07 

WRITE(61*126)GC»OM»TC»PC 

WRITE(61»128)VC 

DO 1 I=1»NH 

1 REA0(60»101) TT( I) »HT( I ) »FMT( n 
DO 2 J»1»NP 

2 READ(60»101 )FMF( J) .CPF( J*l) .CPF( J»2) ♦CPF(J»3» »CPF( J.4) 

DO 3 I«l»NPO 

3 READ (60. 101) FMO ( I ) .CPO( I .1) »CPO( 1.2) »CPO( I .3) »CPO( 1.4) 
REAO(60»100)NO.NT»N.NOR 

WRITE(61.130) 

DO 4 I*1»N0R 

READ (60. 101 ) AO ( I) »G1 ( I ) »0M1( I ) 

0M2( I )=0M1( I ) 

G2( I )aGl( I) 

4 WRITEC61.124) I »A0( I) 

WRITE(61»108) 

DO 5 Isl.N 

READ(60.101) A( I ) .V0{ I ) .Rl( I ) »P1( I ) »AM1( I ) 

READ (60. 10 1 ) AL( I ) .HR( I) »CCF( I ) »AL1( I ) 

R2( I )=R1 ( I ) 

P2( I )=P1( I ) 

AL2< I )*AL1( I ) 

AM2(I)=AM1(I) 

WRITE (61 *125) I » A( I ) *V0( I ) *R1( I ) »P1( I ) 

5 WRITE(61.127) AL ( I ) .HR ( I ) .CCF ( I ) 

AM2(N+1)».0 

AMl(N+l)a.O 

WRTTE(61.108) 

RETURN 

124 F0RMAT(12H ORIFICE NR.I4.6H AREA E15.8.5H (M2)) 

130 F0RMAT(31H ORIFICE DATA FIRST COMPARTMENT) 

101 F0RMAT(5E15.8 ) 
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100 FORMAT(9I3) 

128 FORMATdSH CR I T « V I SCOS I TY *E 13 .6 »6H <G/MS) 

126 FORMATdlH GAS CONST . *£13 .6 »5 h M/0K»7H MOL W.»E13.6»9h CR.TEMP.»E1 
13.6*3H OK.9H CR.PRES. »E13.6 lAM ATM) 

127 FORMAT(22H LEAK AREA»E1S.8(3H M2dOH HYDR .RAO* E13.6 »2H 

1 M.16H RAD.CLEAR/FL.L. »E13.6> 

125 FORMAT(9H COMP. NR* » I 3 dOH OR I F* AREA « E15 * 8 * 3H M2»5H VOL* *E15.8»3H M 
13«8H DENSITY*E15.8*8H KGS2/M4*7H PRES* E15.8.6H KG/M2) 

108 FORMAT dHO) 

END 
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SUBROUTINE OuTPUT(INT) 

C VENTING OF COMPARTMENTS IN SERIES 
COMMON MPRNT »NOR»NPO 

COMMON NH»N.NT*NP.GC»OM»TC»PC»G*TT(40) »HTUO) »FMT(AO) »FMF(AO) *W»NO 
COMMON CPF(40»5) * A ( 5 ) *HR ( 5 ) .PL ( 1 ) *VO ( 5 ) . AL ( 5 ) »CCF<5) »PA1(5) »PZ(1) 
COMMON Pl(5) .Rl(5) .AMKS) .T1.CMK5) ,AL1(5) »G1 ( 5) .OMl ( 5 ) 

COMMON P2(5) »R2(5) .AM2(i>) .T2»CM2(5) »AL2(5) »G2 ( 5 ) .OM2 ( 5 ) 

COMMON PG(5) »CPO(40»5) *FM0(40) »AO(5) tPA » TA »RA * CA » ALT » FM . WM ( 40 . 5 ) 
COMMON TP(40) »OP(40) »WM(40»5) 

WKIT£<61.5C)T1 

50 FORMATdlH TIME E12.5) 

WRI TE (61 .60)ALT.FM.PA.TA.RA 

60 FORMATdlH ALTITUDE ^E12.5»11H MACH NR* E12.5.11H AMB.PR* E12*5 
l.llH AMB.TEMP E12.5.11H AMB*OENS E12*S) 

DPN«Pld)-PA 

1*1 

WRITE(61»61 ) I .PI d > .Rid) .OPN.AMld) 

61 FORMAT (12H COMPARTM. I3.19H PRESSURE E12.5.11H DENSITY 

1 E12.5.11H PRES.DIFF.E12.5.11H MASS FLOW E12.5) 

WRITE(61.65)PAld) .ALld) 

65 F0RMATO4H LEAK PR. E12.5.11H LEAK M.FL.El 

12.5) 

DO 8 I-l.NOR 
PRR = Pld)/PG( I ) 

8 WKIT£(61.66) 1 tPGd ) .Gld ) .OMK I ) .PRR 

66 FORMATdSH ORIF.I3.16H OR I F. PRES. E12 . 5 . 1 IH MASS FLOW 

1 E12.5.11H ORIF.M.NR.E12.5.11H PRES.RAT IOE12.5 ) 

DO 1 1=2. N 
OPN*Pld )-PA 

WRITE(61.63) I .PI (I ) .Ri (I ) .CMld ) .AMI d ) 

1 WRITE (61 .64)PA1 ( I ) .ALl d ) .DPN 

63 FORMAT (12H COMPARTM. I3.19H PRESSURE E12. 5. IIH DENSITY 

1 E12.5.11H MACH NR E12.5.11H MASS FLOW E12.5) 

64 F0RMAT(34H LEAK PR. E12.5.11H LEAK M.FL.El 

12.5.11H PRES.DIFF.E12.5) 

IF( INT-NT)4.6.6 

6 N01=NT+1 
NT1=120 

WRITE(62.80)N01.NT1.N.NOR 
DO 9 1=1. NOR 

9 WRITE(62.81 )AO( I ) *G1 ( I ) .OMK I ) 

DO 2 1=1. N 

WRITE(62.81 )A ( I ) .VO( I ) .Rl ( I ) .Pl( I ) .AMK I ) 

2 WRITE(62.81 )AL( I ) .HRd ) .CCF(I ) »AL1( I ) 

4 DO 3 1=1. N 
AM2 ( I )=AM1 ( I ) 

IF( INT-D5.7.5 

7 AM2( I )=AM1( I )+1.0E-06 

5 P2 ( I )=P1 ( I ) 

R2( I )=R1( I ) 
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CM2U )»CMU I ) 

3 AL2(I)>ALUn 
AM2(N+U»0.0 
AM1(N+1)=0.0 
PH i)».995*pi a ) 
DO 10 l=l*NOR 
G2( I )a61(I ) 

10 0M2(n-0Ml(I) 

RETURN 

80 FORMAT<9I3» 

81 FORMAT(3E15.8) 
END 
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c 


SUBROUTIN E APPT ( K » PX ♦ Z2 .Z3 » A4 ,D1 

cALcuLAtEi If Pl6w li SON I C Oft subsonic 

COMMON MPRNT»NOR»NPO 


COMMON NH»N»NTiNP»GC»OM»TC»PC»G 
AlsPX/ZS 

CALL COEriAl.AIO 
IF.(A1-D)3»3*4 

4 A2»AK*SORT ( 2 .*Z3*Z2 ) *Z4 
A2-A2*(A1*#(1./G) ) 

A3»A1*<M (G-1. )/G) 

IFIA3-1. )5»6.6 


5 A4»A2*S0RT (G* ( 1 .-A3 > / ( G-1 . )) 
GO TO 20 


IN COMPARTMENT APPT* 


6 A4*-A2*SQRT < G*ABS < 1 .-A3 ) / < G-1. )) 

GO TO 20 

3 A2=AK*Z4#SQRT (2.*Z3*Z2) 

A4=A2*( (2./(G+l. ) )**( l./(G-l. ) ) )*SQRT (G/(6+l. J > 
20 RETURN 
END 
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c mi-tmummum c .« or.fke. 

C Zl PRESSURE RATIO. EXTERN. TO INTERN. PRESSURE 

C A INTERN. PRESSURE 

C B INTERN. DENSITY 

C C CLOSE CLEARENCE ORIFICE AREA 

C G SPEC. HEAT RATIO 

C D DISTING. VALUE FOR CHOKING 

C E MYORAUL. RADIUS 

C PC CRITICAL PRESSURE OF GAS 

C TC CRITICAL TEMP. OF GAS 

C R GAS CONSTANT FOR PART. GAS. 

C OM MOL WEIGHT OF GAS 

C AA MASS FLOW THROUGH CLOSE-CLEAR ORIF. 

COMMON MPRNT.NOR.NPO 

COMMON NH»NtNT»NP»RG.OM»TC»PC»G 

VC«7.7*SQRT (OM)*(TC**(-l./6. ) ) »(PC**(2./3. ) )*1.0E-07 
CALL COEF(Zl»AK) 

IF(Z1-D)2»2.1 

1 A2*AK*SQRT J2.*A#8)*C 
A2=A2*(Z1**( l./G) ) 

A3=Z1**( (G-1. )/G) 

IF(A3-1. )3»4.4 

3 A4=A2*SQRT (G*( 1.-A3)/(G-1. ) ) 

L«-l 

GO TO 5 

4 A4=-A2»SQRT (G*ABS (1 .-A3 ) / ( G-1. ) ) 

L=-l 

GO TO 5 

2 A2=AK*C*SQRT (2.*A*B) 

A4*A2»((2./(G+1. ))*»(!. /(G-1. )))*SORT (G/(G+1.)) 

L=1 

5 T=A/(B*RG*9. 78035) 

CALL VISC(T»TC»VC»AMU) 

RE=ABS ( A4»E/ (C*AMU) ) 

CALL CLEAR(RE»RCF*Y) 

A4*A4*Y 

CALL COEF(Zl.X) 

A4=A4*X 
RETURN 
END * 
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SUBROUTINE ATMOS ( HG.P ,T »R ,C1 

AM»28*9&44 

RS=8A7.8118 

A=6378178. 

GO>9.8066S 

GA»1,4 

IF(MG-11000. ) 1»2»2 
1 P8-10332.3 

TP— 6*5/1000. 

TM=288.15 

HB»0. 

GO TO 10 

2 IF(HG-20000. )3»4»A 

3 PB=2307.83 

TP = 0. 

TM=216.65 

HB»11000. 

GO TO 12 

A IF(HG-32000. )5*6»6 
5 PB*568.283 

TP*1./1000. 

TM=216*65 

HB«20000. 

GO TO 10 

6 IF(HG-47000. )7.»8»8 

7 PB=88.513 

TP=2. 8/1000. 

TM=228.65 

HB=32000. 

GO TO 10 

8 IF(HG-52000. )9»13.13 

9 PB*11.309 

TM=270.65 
HB=A7000« 

TP=0. 

GO TO 12 
13 PB=6.016A 

TM=270.65 
HB=82000. 

TP=-2./1000. 

10 A=ALOG(TM/( TM+TP*(HG-HB) ) ) 
P=PB*EXP ( AM*A/ (RS*TP) ) 

GO TO 11 

12 P=PB*EXP (-AM*(HG-HB)/(RS*TM) ) 

11 T=TM+TP*(HG-HB) 
r=AM«P/(RS*T*GO) 

C=SQRT (GA»RS*GO*T/AM) 

RETURN 

END 


SUBROUTINE INTER(N*M>NN»MN»X»Y,XS,YS ) 

DIMENSION X(S0) »Y(50) >XS(14U) »YS(140)*AL(S) 

DO I I»NN.MN 

00 2 L»N.M 

IF«XS( I »-X(L) )3»2»2 

2 CONTINUE 
U«M 

3 L1»L 

IF(L1-N-X)4.4»S 

4 L1»N 

60 TO 8 

5 IF(L1-M+1>6»6»7 

6 Ll=Ll-2 
GO TO 8 

7 Ll*M-3 

8 L2»Ll+3 

00 9 J»L1.L2 
I I»J-L1+1 
AL( I I»=l. 

DO 9 LP=L1»L2 
I J*LP-L1+1 
IF(II-IJ)11»9»11 

11 AL( II )«AL( I I )*( XS( I )-X(LP) )/ (X( J)-X(LP) » 

9 CONTINUE 
SUM =0. 

00 10 IP=Ll»L2 
I I=IP-L1+1 

10 SUM=SUM+AL( I I )*Y( IP) 

1 YS(I)=SUM 
RETURN 
END 


.S uB Ray J IIilS Y lSCJT >TC tY(;:yA> 

TR=T/TC 

IF(TR-1.5)l*lt2 

1 

AN=.9882 
GO TO 9 

2 IF(TR-2. )3>3f4 

3 Al*.462 
ANo.6204 
GO TO 9 

4 IF(TR-3. )5*5»6 

5 Al=.5 
AN=,6480 
GO TO 9 

6 IF ( TR-5* ) 7,7^8 

7 Al=.564 
AN=,5185 
GO TO 9 

8 Al=.631 
AN=#5099 

9 A=( A1*TR)**AN 
A=A*VC 
RETURN 

END 
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SUBROUTI NE CLEAR ( HE »KCF t Y ) 

IF(RE-100.) lfl»2 

1 A>.02 
AN*. 1867 
GO TO 5 

2 IF(RE-500. >3.3.4 

3 A».52 
AN*. 0782 
GO TO 3 

4 A»1.0 
AN»0. 

5 Y»(A+AN*AL0G( l.+RE) )*(1.-EXP (-4.788*SORT (RCF>>) 
RETURN 

END 


S V9RQIJTINE C9£F(A1*AK,> 

CON=1.0 

AB=A1 

IF(AB-1.)1.1.2 
2 AB=1./A1 

CON=-1.0 
1 A=EXP (AB) 

B=l./A 

C*EXP (-3.5*AB) 

AK=( l.-.155*C-.5186*(A-B)/ (A+B) )*CON 

RETURN 

END 


5^UBR0UTINF VFNTCtX»AlC»Y) 

S«4, 5714285 

U *.0101 

W=X**S 

0=W/( ( l.-X)*U+W) 

Y=AK*Q 

RETURN 

END 



PROGRAM VENT 

C VENTING THROUGH A PIPE* NO HEAT ADDITION 
C 

C NH » NUMBER OF CAROS FOR TRAJECTORY DATA 
C N -NUMBER OF COMPARTMENTS ♦ ALSO NUMBER OF LEAKS 
C ONE LEAK PER COMPARTMENT 

C PIPE IS CONSIDERED AS COMPARTMENT NO. 1 

C HT » ALTITUDE AT TIME TT IN METERS 
C FMT • MACH NUMBER AT TIME TT 
C A • ORIFICE AREA BETWEEN N COMP. IN METERS 

C VO « VOLUME OF COMPARTMENTS IN METERS3 

C CPF- PRESSURE COEFFICI.ENT AT MACH N. 

C A » AREA OF NV ORIFICES OF COMP. 1 IN METERS 
C P « PRESSURE IN COMPARTMENT IN (KG/M2) 

C CMl- COMPARTMENT MACH NUMBER 

C ALl* LEAK MASS FLOW 

C R « DENSITY OF N COMP. IN KG.SEC./M4 

C AL -LEAK AREA OF COMP.N IN METERS 
C HR -HYDRAULIC RADIUS OR MEAN RADIUS IN MET&R 

C TC » CRITICAL TEMPERATURE OF GAS 

C GC - GAS CONSTANT OF GAS (M/DEG. K) 

C PC - critical pressure of gas in ATM(KG/CM2) 

C OM » MOL WEIGHT OF PARTICULAR GAS 
C CCF- RADIAL CLEARANCE OVER FLOW LENGTH 

C WI -ADDITIONAL LOSS COEFFICIENT FOR PIPE* IF NONE ENTER 0. 

C IN PROPER FORMAT 

C NUMBER 1 COMP. IS THE COMP. THAT VENTS INTO THE ATMOSPHERE 

COMMON MPRNT 

COMMON NH.N»NT*NP*GC.OM»TC.PC*G.TT<40) *HT(40) * FMT (40) »FMF(40) »W»NO 
COMMON CPF (40. 5 ) » A ( 5 ) . HR ( 5 ) . PL ( 1 ) * VO ( 5 ) * AL ( 5 ) *CCF( 5 ) * PAl ( 5 ) . Pi ( 1 ) 
COMMON P1(5).R1(5) .AMK5) *T1.CM1(5) *AL1(5) 

COMMON P2( 5) »R2(5) »AM2(5) .T2.CM2(5) *AL2(5) 

COMMON WH(40»5) *TP(40) *OP(40) »EP * ER * EM *EMX *PA * TA.RA.CA » ALT » FM 
COMMON WM(40*5) 

DO 5 1-1*5 
AMI ( I >»0.0 

5 AM2( I )»1.0E-07 
CALL INPUT 

D-(2./(G+l. n**(G/(G-l. ) ) 

DO 13 I-NO.NT 
1 1=1-1 
Tl-Il 
T2-T1-1.0 

CALL TRAJ(ALT»FM»HG*PA*TA*RA*CA) 

IF( 1-1)4.11.4 
4 DO 15 L-1.40 
IF(L-1)6.6.7 

6 DO 8 K-l.N 

8 AM1(K)=AM2(K) 

7 DO 20 K-2.N 
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IF(K-2)1»1.2 

1 CONTINUE 

CALL PARTB( K*L# I ) 

GO TO 23 

2 CONTINUE 

CALL PARTA( 1 «K»L*0) 

23 WH(L»K)«P1(K-1> 

20 WM(L.K)»AMl(IO 
IF(L-1 ) I0.10t9 

9 DO 3 <«1»N 

IF(MPRNT>18il8.19 

19 WRITE(6l«91)K»WM<L»K) .WM(L-1*K) 

18 IF(ABS(WM(L*K)-WM(L-l.ia ) -1 .OE-07 ) 3 1 3 » 10 

3 CONTINUE 
GO TO 12 

10- IF(MPRNT)21.21»22 

22 WRITE(61»90)L*L*AM1( 1) • AMI ( 2 ) *AM1 ( 3 ) »AM1 ( 4 ) 

21 IF(L-3) 15*16 1 16 

16 DO 17 K»3»N 

17 PKK-l )«WH< L»K)-<WM(L*K)-WM(L-1*K) ) * ( WH ( L *»0-WH ( L-1 tK ) ) / ( WM ( L *10 
12,.0#WM(L-1.KJ+WM(L-2»K) ) 

15 CONTINUE 
GO TO 12 

11 EM=0.0 
EMX=0.0 

DO 14 IP=1*N 
CMl ( IP) »+0*0 
14 ALl(IP)*+0.0 

12 CALL OUTPUT! I ) 

WRITE(61»555) 

13 CONTINUE 
STOP 

555 FORMAT (IHO) 

90 F0RMAT(2I4*8E15.8) 

91 FORMAT! I4*5E15. 8) 

END 


SUBROUTINE PART A ( I ♦ K » L *0 > 

DIMENSION T (^1 ) .0(41 ) 

COMMON MPRNT 

COMMON NH.N»NT.NP»GC.0M.TC»PC»G.TT(40> »HT(40) .FMT(40) » FMF ( 40 ) . W .NO 
COMMON CPF (40. 5) . A ( b ) .HR ( 5 ) .PL ( 1 ) . VO ( 5 ). ♦ AL < 5 ) .CCF ( 6 ) .PAl ( 1> ) »PZ ( 1 ) 
COMMON PK5) .RK5) .AMl(ij) .T1.CMK5) .AL1(1>) 

COMMON P2(5) .R2(5) .AM2(5) .T2»CM2(5) .AL2(5*) 

COMMON WH(40.5) .TP(40) .OP(40) .EP.ER.EM.tMX.PA.TA.RA'iCA.ALT.FM 
T( l)«(P2(IO+Pl(K) )/2.0 

DO 30 M=1.40 

A1=(P1(K-1)+P2(K-1) )/2.0 
R1(K)=R2(K)*( ( T(M)/P2(ia )**(!. /G) ) 
call APPT(K.A1.R1(<) .T(M) .a (IO.A4.D) 

IF(AL(K) >27.27.33 
27 AL1(K)=0.0 

GO TO 28 
33 Z1=PA/T(M) 

CALL LEAK(Z1.T(M) . R1 ( K ) . AL ( K) .0 .HR ( K) .CCF(lO .ALKK) ) 

28 • Sl=R2(K)*VO(K) 

S2= (2.0*( T(M) /P2(K) ) -1 . 0 ) ** < 1 . /G ) 

S3=S1*( 1.0-S2 )+AMl (K1 )-ALl (K) 

CMX2=2.0*( ( T (M) /A1 )**( (G-1. ) /G)-1.0 ) / (6-1. > 

CM1(K)=S0RT(ABS(CMX2) ) 

IF (CMl (10-1.0)6. 6.7 
7 CM1(K)«1.0 

6 0(M)=A4-S3 

IF(M-l) 54.54.35 
54 T ( M+1 ) =T (M) *1 .005 
GO TO 62 

35 T(M+1)=T(M)-0(M)*(T(M)-T(M-1) ) / ( 0(M ) -0(M-1 ) ) 

62 P1(K)=2.0»T(M+1)-P2(<) 

AMI (K)=S3 
IF(MPRNT-2) 1.1.3 

3 WRITC(61.90)M.T(M).A1.A4.A(IO.AL1(K) .S2 .S3 . T (M+1 ) 

90 F0RMAT( I3.8E15.8) 

1 IF(ABS (0(M) )-1.0E-07)22.30.30 
30 CONTINUE 

WRITE(61.99) 

22 IF(MPRNT)4,4.5 

5 WRITE(61.80)I.K.L.M.D.T(M).P1(IO.S3.0(M).R1(K) 

4 RETURN 

30 F0RMAT(4I4.6E15.8 ) 

99 FORMAT (17H PARTA TOLERANCE) 

END 


■4 
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SUBROUTINE PARTB ( K1 . LI » I ) 

DIMENSION PM(41) 

COMMON MPRNT 

COMMON NH*N*NT*NP*GC.OM»TC#PC»G#TT(40) »HT(40) ♦FMTUO) »FMF(40) *W»NO 
COMMON CPF(40*5 ) . A ( 5 ) tHR < 5 ) ,PL ( 1 ) ♦ VO ( & ) i AL ( 5 ) »CCF ( 5 ) »PA1 ( 5 ) ♦PZU ) 
COMMON PK5) »R1(5) *AM1(5) .T1.CMK5) *AL1(5) 

COMMON P2«5) iR2(5) .AM2(5) iT2fCM2(5) tAL2(5) 

COMMON WH(40.5) .TP(40) tOPUO) . EP i ER * EM . EMX * PA . TA . RA *C A » ALT ♦ FM 
COMMON WM(40*5) 

K»K1-1 
NU = 0 
NS»0 

DO 1 L=l»40 
IF(MPRNT-1 ) 2»2f 10 

10 WRITE(61.555) 

WRITE(61»100)L*K 

2 IF(L-2)6.3.5 

6 P.M(L)aPKKl) 

GO TO 5 

3 PM( L) =.998»PM(L-1 ) 

5 P1(K1)=PM(L) 

TL=P1 (K1 ) /( R1 (K1 )*GC*9t8) 

CALL ENTR(K.SV) 

PMl L) =P1 ( K1 ) 

AMKK1)=AM1(K) 

HPO=EM 

TL»TL/« l.+(G-l. )*EM *EM /2.) 

call FRICTIAMinC ) . A (,K, ) i T L t FR . V I S *RN ) 

FF»G*(FR*PL(K)/ (2.*HR(K) ) +W ) 

CALL PIPE(FF,CE»SV»HP.GLIM) 

IF(MPRNT-l) 11.11.12 

12 WRITE(61.100)L»K.L1.P1(K) » EP ♦ EM. EMX . PAl ( K > .AMKIO .AMIUCI) .HP 
WRITE(61.100)L.K1.L1.P1(K1).R1 (k;).P2(K1).P2(K) 

11 DE=Pl(K)-PAl(lO 
DEM=HPO-HP 

CALL ITER(NU.NS.PM(L) .DE . DEM . PM ( L+1 J .MS8.MS0) 

IF(MPRNT-l) 13.13.14 

14 WRlTE(61.i01)L.K.AMl(Kl) .PM( L ) .PM( L+1 ) .DE.DEM 

13 IF(EMX-1. )17.7.7 

17 IF(MSB)4.8.4 

4 IF«ABS(DE)-(1.0E-08*P1(K1) > )8.8.1 

7 IF(DE+(2.0E-08*P1(K1) ) )1.9.9 

9 IF(ABS(OEM)-1.0E-07)8.8.18 

18 IF(MSO)1.8.1 
1 CONTINUE 

8 WM(L1.K)=AM1<K) 

WM( LI .XI ) =AM1 (XI ) 

CM1(X1)=EM 

CMl (X 1 =EMX 
IF(MPRNT) 15.15.16 
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16 WRITE(61*102)L1»L 

15 RETURN 

101 FORMAT(2I4.6E15.8 

100 FORMAT(3I4»8E15.8 

102 F0RMAT(4I4.8£15.8 
555 FORMAT (IHO) 

END 


*MSB*MS0*AM1(K) »P1 ( Kl ) »P2 ( K1 ) 

) 

) 

) 
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SUBROUTINE INPUT 
COMMON MPRNT 

COMMON NH»N»NT.NP»GC»OM»TC.PC.G.TT(40) »HT(40) »FMT(40) ♦ FMF ( 40 » ♦ W #N0 
COMMON CPF(40»5) »A(5) »HR(5) »PL(1 ) ♦ VO ( 5 ) » AL { 5 ) »CCF ( 5 ) »PA1 ( 5 ) »PZa) 
COMMON P1(5)»R1<5) »T1»CM1(5) *ALl(5) 

COMMON P2(5) »R2(5) »AM2(5) »T2»CM2(5) »AL2(5> 

COMMON WH(40»5) »TP(40) .0P(40) » EP »ER * EM» EMX »PA ♦ TA»RA»CA» ALT ♦ FM 

READ(60*100)NM»NP»MPKNT 

READ(60»101 )OM»TC.PCtOC»C 

VC»7.7*SQRT (0M)*(TC**(-l»/6. ) )*{PC**(2*/3. ) )*1.0E-07 

WRITE(61»126)GC»0M»TC»PC 

WRITE(61*128»VC 

DO 1 I*1*NH 

1 READ(60»101) TT( n tHTm »FMT( I ) 

00 2 J«1»NP 

2 READ! 60* 101 ) FMF (J) .CPF ( J . 1 ) , CPF ( J» 2 ) »CPF(J»3) »CPF(J*4) 
READ(60»100)NO.NT*N 

WRITE(61*130) 

READ (60 *101 ) A( 1 ) »HR( 1) »PL( 1) »R1 (1) »P1(1) 

WRITE(61*124) A( 1 ) *HR< 1 J *PL( 1) 

WRITE{61»108) 

DO 5 I»2*N 

READ (60* 10 1 > A( I) *V0( I ) *R1( I ) *P1( I ) *AM1( I) 

REA0(60»101) AL( I ) *HR( I) *CCF( I ) *AL1( I ) 

R2(n=Rl(I) 

P2( I )aPl( I ) 

AL2( I )«AL1( I ) 

AM2( I )=AM1( I ) 

WRITE(61*125) I* A( I ) »V0( I) *R1( I) *P1( I ) 

5 WRITE(61*127) AL ( I ) *HR ( I ) »CCF ( I ) 

AM2(1)=AM1(2I 
EP=P1 (2) 

ER=R1 (2) 

READ(60*101 )W 
WRITE(61*108) 

RETURN 

124 FORMATdOH PIPE AREA » E15 . 8 . lOH HYDR.RAD. *E1S.8 .8H PIPE L..E15.8) 
130 FORMAT! 13H PIPE SECTION) 

101 FORMAT! 5E15. 8) 

100 F0RMAT(9I3) 

128 F0RMAT(15H CR I T . V I SCOS I TY *E 13 .6 * 6H KG/MS) 

126 FORMATdlH GAS CONST. *E13.6*5H M/0K.7H MOL W.*E13.6*9H CR.TEMP..E1 
13.6.3H 0K»9H CR.PRES. *E13 .6 *4H ATM) 

127 FORMAT(22H LEAK AREA *E1 5 . 8 » 3H M2.10H HYDR .rRAD. E 13. 6 » 2H 

1 M.16H RAD.CLEAR/FL.L. *E13.6) 

126 F0RMAT(9H COMP. NR. » I 3 * lOH ORIF.AREA*E15.8*3H M2.5K VOL. *E15 . 8 * 3H M 
13*8H DENSITY.E15.8.8H KGS2/M4.7H PRES. E15.8.6H KG/M2 ) 

108 FORMAT! IHO) 

END 
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50 

60 


61 


62 


1 

63 

6 ^ 


6 


2 

4 


7 

5 


3 
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SUBROUTINE OUTPUT(INT) 

COMMON MPRNT 

COMMON NH*N»NT.NP»GC»OM*TC»PC»G.TT(40) *HT(40) »FMT(40) »FMF(40) »W*N0 
COMMbN CPF(40»5) »A(5) fHR ( 5 ) .PL ( 1 ) ♦ VO ( 5 ) » AL ( 5 ) »CCF I 5 ) »PA1 ( 5 ) »PZ ( 1 ) 
COMMON P1(5)»R1(5) .AMI ( 5 ) ♦ T 1 »CM1 ( 5 ) .ALK5) 

COMMON P2(5).R2(5) .AM2 ( 5 ) . T2 .CM2 ( 5 ) .AL2(5) 


COMMON WH(40.5) .TP(40) .OP(40) .EP.ER.EM.EMX.PA»TA.RA.CA.ALT.FM 
COMMON WM(40»5) 

WRITE(61.50)T1 
FORMATdlH TIME E12.5) 

WRITE(61.60) ALT.FM.PA.TA.RA 

FORMATdlH ALTITUDE E12.5.11H MACH NR. E12.5.11H AMB.PR. E12.5 
l.llH AMB.TEMP E12.5.11H AMB.DENS E12.5) 
■wRITE<61.61)EP.ER.EM.AMld) 

F0RMATI34H COMPARTM. PIPE ENTR.PR. E12.5.11H ENTR.DENS.El 

12.5.11H ENTR.M.NR.E12.5.11H ENTR.M.FL.E12.5) 

WRITE(61.62)Pld) .Rim .EMX.AMld) 

FORMAT(34H EXIT PR. E12.5.11H EXIT. DENS El 

12.5.11H EXIT M.NR.E12.5.11H EXIT M.FL.E12.5) 

DO 1 1=2. N 
DPN = P1 ( n-PA 

WRITE (61.63)1. Pim. Rim .CMl ( I ) .AMI ( I ) 

WRITE(61.64)PA1( I) .ALld ) »OPN 


FORMAT(12H COMPARTM. I3.19H PRESSURE E12.5.11H DENSITY 

1 E12.5.11H MACH NR E12.5.11H MASS FLOW E12.5) 

F0RMAT(34H LEAK PR. E12.5.11H LEAK M.FL.El 

12.5.11H PRES.OIFF.E12.5) 

IF( INT-NT)4.6.6 

N01=NT+1 

NT1=120 

WRITE(62.80)N01 .NTl.N 

WRITE(62.81 )Ad) .HRd) .PLd) .R1 d) .Pld) 

DO 2 1=2. N 

WRITE(62.81 )A(I > .VOd ) .R1 ( I ) .Pl( I ) .AMK I ) 

WRITE(62.81 )AL(I ) .HRd) .CCF( I) »AL1( I ) 

WRITE(62.81)W 
DO 3 1=1. N 
AM2( I )=AM1( I ) 

IF( INT-D5.7.5 

AM2 < I ) =AM1 ( I ) +1 .OE-06 

P2d )=P1( I ) 

R2( I )=R1( I ) 

CM2( I )=CM1( I ) 

AL2( I )=AL1( I ) 

AM2(N+1)=0.0 
AM1(N+1)=0.0 
EP=.99*EP 
RETURN 
FORMAT! 91 3) 

FORMAT(5E15.8) 

END 
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SUBROUTINE APPT ( fe ♦ PX « Z2 » Z3 « > A4 « D » 

C CALCULATES IF FLOW IS SONIC OR SUBSONIC IN COMPARTMENT APPT. 

COMMON MPRNT 

COMMON NH»N»NTtNP»GC»OM»TC»PC»G 

A1=PX/Z3 

CALL COEF(Al.AK) 

IF(A1-0)3*3»A 

4 A2»AK*SQRT ( 2 . »Z3*22 I *Z4 
A2«A2*(A1**( l./G> > 

A3*A1**( (G-1. ) /G) 

IF(A3-1.)5.6.6 

5 A4*A2»SQRT ( G* ( 1 .-A3 ) / ( G-1 . ) ) 

GO TO 20 

6 A4«-A2*SQRT ( G*ABS (1 .-A3 ) / ( G-1. ) ) 

GO TO 20 

3 A2=AK*Z4*SQRT (2.*Z3*Z2) 

A4=A2»( (2./ (G+1. ) )**( 1./ (G-1. ) > )*SQRT (G/(G+1. ) ) 

20 IF(A1-1. ) 1.1»2 

2 IF(A4) 1.1.8 
8 A4=-A4 

1 RETURN 
END 



SUBROUTINE LEA K(Z1»A>B»C.D.E.RCF>A4) 

c calculates if flow is sonic or subsonic in orifice. 

C 21 PRESSURE RATIO. EXTERN. TO I/^TERN. PRESSURE 

C A INTERN. PRESSURE 

C 8 INTERN. DENSI TY 

C C CLOSE CLEARENCE ORIFICE AREA 

C G SPEC. HEAT RATIO 

C D DISTING. VALUE FOR CHOKING 

C E HYORAUL. RADIUS 

C PC CRITICAL PRESSURE OF GAS 

C TC CRITICAL TEMP. OF GAS 

C R GAS CONSTANT FOR PART. GAS 

C OM MOL weight OF GAS 

C A4 MASS FLOW THROUGH CLOSE-CLEAR ORIF. 

COMMON MPRNT 

COMMON NH.N.NT.NP.RG.OM.TC.PC.G 

VC=7.7*SORT (OM)*(TC**l-l./6. ) ) * ( PC** < 2 . /3 . ) )*1.0E-07 
CALL COEF(Zl.AK) 

IF(21-D)2.2»1 

1 A2=AK*SORT (2.*A*8)*C 
A2=A2*(Z1**(1./G) ) 

A3=Z1**( (G-1. ) /G) 

IF(A3-1. >3,4,4 

3 A4=A2*SQRT «G*( 1.-A3)/(G-1. ) ) 

L=-l 

• GO TO 5 

4 A4=-A2*SQRT ( G*ABS ( 1 .-A3 ) / ( G-1 . ) ) 

L = -l 

GO TO S 

2 A2=AK*C*SQRT (2.*A*B) 

A4 = A2*( ( 2./ (G+1. ) > **( 1./ ( g-1. ) ) )*SQRT (G/ (G+1. ) ) 

L=1 

5 T=A/(B*RG*9.7803S) 
call VISC(T,TC,VC,AMU) 

RE=ABS (A4*£/(C*AMU) ) 

CALL CLEAR! RE, RCF.Y) 

A4=A4*Y 

CALL COEF(Zl,X) 

A4=A4*X 

RETURN 

END 




76 



SUBROUTINE P I PE ( FF » VX » EV » H .GL IM > 

C H = ENTRANCE MACH NUMBER FOR CIhOKED FLOW EMX-1 
C EM =MACH NUMBER AT ENTRANCE OF PIPE. 

C EMX=MACH NUMBER AT EXIT OF PIPE 
C G = RATIO OF SPEC. HEAT 

C FF = PIPE FRICTION FACTOR 

C EP = PRESSURE AT ENTRANCE 

C £R * DENSITY AT ENTRANCE 

C PX -PRESSURg AT EXIT 

C VX » VELOCITY OF SOUND AT EXIT 

C £V * VELOCITY OF SOUND AT ENTRANCE 

C R1 = DENSITY AT EXIT 

DIMENSION 0(41) »X(41) 

COMMON MPRNT 

COMMON NH»N»NT.NP»GC.OM*TC»PC»G.TT(40) »HT(40) »FMT(40) »FMF(40) »W»NO 
COMMON CPF(40»5) . A ( 5 ) .HR ( 5 ) »PL ( 1 ) » VO ( 5 ) » AL ( 5 ) »CCF(5) »PA1(5) »PZ(1) 
COMMON PK5) .Rl<5) »AM1(5) .T1.CMK5) »AL1(5) 

COMMON P2(5) »R2(5) »AM2(5) »T2.CM2(5) »AL2(5) 

COMMON WH(40.5) »TP(40) »OP(40) . EP .ER » EM» EMX 
I J=1 

G1=(G+1. )/2. 

G2=(G-l.)/2. 

B=EM#*2 

A2«FF-( l./B)-Gl*AL0G(ABS(8/( l.+G2*B) ) ) 

A1=-1.+G1»AL0G(G1) 

IF(A2-A1)8.13»13 
13 EMX=1. 

I J = 2 

8 0(1)=1. 

GO TO 10 
35 0( 1 )=X1 

10 DO 34 1=1.40 
GO T0( 1.2) . I J 

1 Z1=AL0G(B»( l.+G2*D( I ) )/(D( I )*(1.+G2*B) ) ) 

X1=1./(1./B-FF+G1*Z1) 

X(I ) =X1-D(I ) 

EMX=SQRT ( ABS(D( I ) ) ) 

F=EMX*«2 
GO TO 11 

2 21 = AL0G(ABS(0( I )*G1/ ( l.+G2*0( I ) ) ) ) 

Xl=( l./( 1.+FF-G1*Z1) ) 

X( I )=X1-D( I ) 

H=SQRT(ABS(D( I ) ) ) 

11 IF< 1-1 ) 15.15.7 

15 D(2)=.l 

GO TO 34 

7 0( I+1)=0( I )-X( I )-»(D(I )-D( I-l) ) / (X( I )-X( I-l) ) 

IF(ABS(X( I ) )-1.0E-08)6.6.34 
34 CONTINUE 
GO TO 35 



6 GO TO(21.14) flJ 
21 IJ»2 

GO TO 8 

14 IF( EMX-1. ) 17.3*3 
17 XX»1.+G2*F 

XY=1.+G2*B 

PI ( 1 ) =EP*SORT ( ABS(B*XY/ (F*XX) ) ) 

VX»EV*SQRT(ABS(XY/XX) ) 

Rl( 1)=£R»»SQRT<ABS(B*XX/(F*XY) > ) 

2 1-H+H 

Xl=l.+G2*21 

X2=1.+G2*B 

GL I M=G*EP*SQRT ( Z 1*X1 / ( B*X2 ) ) #H*A ( 1 ) / ( EV*SQRT ( X2/X1 ) ) 
GO TO 4 

3 B1=H*H 
XX=1.+G2*B1 
XY=1.+G2*B 

EP=EP/SQRT(B1#XX/(B*XY) ) 

ER=ER*SQRT(B«XX/(B1*XY) ) 

EV=EV*SQRT( XY/XX) 

PI { 1 ) =EP*H*SQRT (XX/Gl ) 

VX = EV*SORT( XX/Gl ) 

R1 C 1 ) =£R*H*SQRT(G1/XX J 
EM=H 

GLIM=G*EP»H*A( 1 ) /EV 

4 RETURN 
END 
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SUBROUTINE I TER ( NU .NS »PM »DE » DEM,PM1 »MSUB ,MSON ) 

DIMENSION 0PX(40) »AMA(40) »OPE ( AO > »AMB (40 ) 

COMMON MPRNT 

COMMON NH.N»NT»NP»GC»OM»TC»PC»G»TT(40) »HT(40) »FMT(40) »FMF(40) »W»NO 
COMMON CPF (40. 5) ♦ A ( 5 ) .HR ( 5 ) »PL ( 1 ) .V0(5) .AL(5) .CCF( 5 ) »PA1 ( 5 ) »Pi( 1) 
COMMON PK5) .Rl(5) .AMI ( 5 » . T 1 .CMl ( 5 ) »AL1(5) 

COMMON P2( 5) .R2 (5) .AM2(5) .T2.CM2(S) »AL2( 5) 

COMMON WH(40.5) .TP (40) .OP (40) .EP.ER.EM.EMX 

MSUQaJL 

MS0N=1 

IF(EMX-1. ) 1.3.3 

1 NU»NU+1 
ANV=NU 
DPX(NU)»0E 
AMA(NU) »PM 
IF(NU-2)6.2.2 

2 AA=A8S(AMA(NU)-AMA(NU-1) ) 

IF(AA-1.0E-07*P1(2) )9.12.12 

12 PM1=AMA(NU)-0PX(NU)*(AMA(NU)-AMA(NU-1) ) / ( OPX ( NU ) -OPX ( NU“1 ) ) 

GO TO 5 

9 MSUB=0 

IF(ABS(DE)-(1.0E-07*P1(2) ) )8»8.15 

15 IF(ABS(0EM)-2.0E-06)16.16.8 

16 MSUB=1 
EMX=1.0 

3 NS=NS+1 
ANS=NS 
OPE(NS)=OEM 
AMB(NS)=PM 
IF(NS-2)8.4.4 

4 AB=ABS( AMB( NS)-AMB(NS-1 ) ) 

IF(AB-1.0E-07#P1(2) ) 13.14.14 

13 MS0N=0 
GO TO 8 

14 PMl =AMB ( NS ) -OPE ( NS ) * ( AMB ( NS ) -AMB ( NS-1 ) ) / ( OPE ( NS ) -OPE ( NS-1 ) ) 
IF(PM1-AM8(NS) ) 10.11.10 

11 PM1=AM8(NS) *( ANS+0.010)/ANS 
10' IF(ABS(0EM)-1.0E-06)7.7.5 

7 IF(DE+1.0E-07)1.5.5 

6 PM1=(PM+P2(2) )/2.0 

GO TO 5 

8 PM1=.9995*PM 

5 RETURN 
END 
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SUBROUTINE ENTRtF>»SV> 

C AMl= MASS FLOW FROM UPPER COMPARTMENT 

C PI « MEAN PRESS. OF UPPER COMPARTMENT "TOTAL PRESSURE 

C Rl " MEAN DENSITY OF UPPER COMP, i TOTAL DENSITY 
C G = SPEC. HEAT RATIO 

C A * PIPE AREA . 

C EP * PRESSURE IN PIPE ENTRANCE (MEAN) 

C EM « MACH NUMB. AT ENTR* OF PIPE 

C ER - DENSITY AT ENTR. OF PIPE 

C SV = VELOCITY OF SOUND IN THE ENTRANCE 
COMMON MPRNT 

COMMON NH»N.NT.NP.GC»0M.TC.PC»G»TT(40) »HT(40) »FMT(40) »FMF ( 40 ) » W »NO 
COMMON CPF(40.5 ) »A ( 5 ) .HR ( 5 ) »PL ( 1 ) . VO ( 5 ) » AL ( 5 ) »CCF( 5 ) »PA1 ( 5 ) »PZ( 1) 
COMMON PK5) .Rl(5) .AM1<5) .T1.CMK5) »AL1<5) 

COMMON P2(5) »R2(5) »AM2(5) .T2.CM2(5) »AL2(5) 

COMMON WH(40.5) »TP(40) »OP(40) .EP.ER.EM.EMX.PA.TA.RA.CA.ALT »FM 
D=(2./(G+1. ) )**(G/(G-1. ) ) 

K1*K+1 
DO 2 1=1.10 
IF(AL(K1) )11.11.12 

11 AL1(K1)=.0 
GO TO 13 

12 Z1=2.0*PA1(K1)/(P2(K1)+P1(K1) ) 

TM=(P2(K1)+P1(K1) )/2.0 

CALL LEAK(21.TM.R1(K1 ) .AL(Kl) .D.HR< K1 ) .CCFCKl ) .ALl (K1 ) ) 

13 S1=R2(K1 )*V0(K1) 

AM1(K)=AM1(K1+1 )-ALl(Kl)+Sl*( l.-( ( PI ( K1 ) /P2 < K1 ) )**( l./G) ) ) 

RKKl )=R2(K1 )*(!.+ (AMI (Kl + 1 ) -ALl (ta )-AMl(K) )/Sl) 

A0=( (AM1(K,)/A(K) )**2>/(G*Pl(Kl)*Rl(Kl) ) 

B=.739S8*G 

A1=(2.*A0»B-1. ) / (G-1.-2.*B*B*A0) 

A2=(2.*A0)/ (G-1.-2.*B*B*A0) 

FMP=A1+SQRT (ABS(A2+A1*A1 ) ) 

RPT=( l.+(G-l. )*FMP/2. )**(G/( l.-G) ) 

PT2=P1 (K1 )/ (RPT*( 1.0+B*FMP) ) 

PTR=1.0/(RPT*( 1.0+B*FMP) ) 

£P=rPT*PT2 

ER=R1 (K1 )*PT2*( ( l.+(G-l. )*FMP/2. )**( 1./ ( l.-G) )) /PI (K1 ) 
EM=SQRT(ABS(FMP) ) 

SV=SQRT(ABS(G*EP/ER) ) 

IF(PTR-1. )6.6.1 

6 ir(AO-( 1.0+G) /(2.0*( ( l.+B)**2) ) )3.3.1 

1 P1(K1) = (P1(K1)+P2(ia) )/2.0 

2 CONTINUE 

3 IF (MPRNT-2) 4.4.5 

5 WRITE (61 ..lOOJPTR.FMP.SV.EP.ER.Al.AO.PKlCl) 

4 RETURN 

100 FORMAT( 8E15.8 ) 

END 
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SUBROUTINE FR 1 CT ( Y » AC » T » FR » A t R > 

DIMENSION F(41) 

COMMON MPRNT 

COMMON NH»N*NT*NPfGC»0M»TC»PC»6 
C Y “MASS FLOW 

C AC “PIPE AREA 

C T “GAS TEMP. 

C TC “CRIT, TEMP. 

C VC -CRIT. Vise. 

C FR “FRICTION FACTOR 

VC“7.7*SORT (0M)*(TC**(“l./6. ) ) * ( PC** ( 2 . /3 . ) )*1.0E-07 
CALL VISC(T»TC*VC.A) 

0-SQRT (A.*AC/3.1A16) 

R“A3S (Y*0*9.8/(AC*A)) 

IF(R)4.4»5 

4 FR“.0 
GO TO 3 

5 F( 1 >«.05 

DO 1 I-1.40 

F( I+l )“!./(( .86858*ALOG(ABS ( R*SQRT ( F ( I ) ) ) ) - . 8 > **2 ) 
FR=F( I+l ) 

IF( I-m.1,2 

2 IFIABS (F(I+1)-F(I) )-1.0E“07)3»l»l 
1 CONI INUE 

IF( .86858*ALOG(ABS(R*SQRT(FRT) ) -.8)6*6 *7 

6 FR“10.0 
GO TO 3 

7 WRITE(61.100) 

3 RETURN 

100 FORMAT (3X.15HFRICT TOLERANCE) 

END 



SU6H0UTINE TRAJ< ALTtFM.HG.PA»TA«RA*CA) 

PA*AMBIENT PRESSURE 
TA»TEMPERATURE 
RA=OENSITY 
CA*SPEED OF SOUND 

PAl (N)»PKESSURE OUTSIDE PIPE AND LEAKS AT TIME T1 
DIMENSION Wl(AO) 

COMMON MPRNT 

COMMON NH»N»NT iNP »GC«0M«TCiPC*6tTT(40) »HT(40) »FMT(40) *FMF(40) »W»N0 
COMMON CPF {40. 5 ) . A ( 5 ) .HR ( 5 ) . PL ( 1 ) ♦ VO ( 5 ) » AL ( 5 ) »CCF { 5 ) . PAl ( 5 ) »P2U) 
COMMON Pl(5) .Rl(5) .AMK5) .Tl.CMKb) .ALK5) 

COMMON P2( 5) .R2(5) .AM2(5) .T2..CM2(5) »AL2(5> 

CALL INTER! 1 .Nl-1 . 1 . 1 . TT .HT . T 1 .ALT ) 

CALL INTER! l.NH.i .1. FT.FMT. ri.FM) 

DO 1 I=1.N 
DO 2 J»1.NP 
WU J)=CPF!J.I ) 

CALL INTER! 1 . NP . 1 . 1 . FMF . W1 .FM.PAl ( I ) ) 

CONTINUE 

HG=ALT<»6.3 78178E+06/ I6.378178E+06+ALT) 

CALL ATMOSIHG.PA.TA.RA.CA) 

00.3 I=1.N 

PAl ! I )=PA+RA*PA1( I )*{ !FM*CA)*»2)/2«0 

RETURN 

END 



PROGRAM VENT 

C VENTING THROUGH A PIPE. HEAT ADDITION TO PIPE 
C *##*##•**##*»#»*##*#*###*«•###**#*»#*#***#****#*********** 

C NH = NUMBER OF CARDS FOR TRAJECTORY DATA 
C N “NUMBER OF COMPARTMENTS. ALSO NUMBER OF LEAKS 
C ONE LEAK PER COMPARTMENT 

C PIPE IS CONSIDERED AS COMPARTMENT NO. 1 

C HT = ALTITUDE AT TIME TT IN METERS 
C FMT “ MACH NUMBER AT TIME TT 
C A » ORIFICE AREA BETWEEN N COMP. IN METERS 

C VO » VOLUME OF COMPARTMENTS IN METERS3 

C CPF“ PRESSURE COEFFICIENT AT MACH N. 

C A “ AREA OF NV ORIFICES OF COMP. 1 IN METERS 
C P = PRESSURE IN COMPARTMENT IN (KG/M2I 

C CM1= COMPARTMENT MACH NUMBER 

C AL1= LEAK MASS FLOW 

C R « DENSITY OF N COMP. IN KG.SEC./MA 

C AL “LEAK AREA OF COMP.N IN METERS 

C HR “HYDRAULIC RADIUS OR MEAN RADIUS IN METER 

C TC “ CRITICAL TEMPERATURE OF GAS 

C GC “ GAS CONSTANT OF GAS (M/OEG.K) 

C PC “ CRITICAL PRESSURE OF GAS IN ATM(KG/CM2) 

C OM “ MOL WEIGHT OF PARTICULAR GAS 

C CCF« RADIAL CLEARANCE OVER FLOW LENGTH 

C WI “ADDITIONAL LOSS COEFFICIENT FOR PIPE. IF NONE ENTER Q. 

C IN PROPER FORMAT 

C NUMBER 1 COMP. IS THE COMP. THAT VENTS INTO THE ATMOSPHERE 

C OT = HEAT ADDITION ALONG TRAJECTORY 

C Q » TOTAL HEAT ADDITION TO THE PIPE AT TIME T1 

COMMON MPRNT 

COMMON NH»N.NT.NP.GC.OM.TC.PC.G.TT(40) »HT( AO) .FMTIAO) .FMF(40) »W.NO 
COMMON CPF(AO.S) . A < 5 ) .HR ( 5 ) .PL ( 1 ) . VO < 5 ) . AL ( 5 ) .CCF ( 5 ) » PAl ( 5 ) .PZ ( 1 ) 
COMMON Pl( 51 .R1 (5) .AMK5) .T1.CMK5) .ALK5) 

COMMON P2(5) .R2(5) .AM2(5) .T 2. CM2 (5) .AL2(5) 

COMMON WH(A0.5) .TP(AO) .OP(40I . EP . ER » EM. EMX .PA . TA .RA .CA . ALT . FM 
COMMON WM(40.5) .QT(AO) .Q.AM3(5) .P3(5) 

DO 5 I “1.5 
AMKI )“0.0 

5 AM2( I)=1.0E-07 
CALL INPUT 

D“(2./(G+1. ) )»*(G/(G-1. ) ) 

DO 13 I=NO.NT 
I 1“I-1 
T1=I1 
T2“T1-1.0 
CALL TRAJ 
IF( 1-1)4.11.4 
4 DO 15 L“1.20 
IF(L-1)6.6.7 

6 DO 8 K=1.N 
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P1(K)*2«0*P2(K)-P3(K) 

8 AMl(IOa2.0#AM2(K)-AM3<K) 

7 DO 20 K=2»N 

IF(K-2) 1»1»2 

1 CONTINUE 

CALL PARTB(K.L»I) 

GO TO 23 

2 CONTINUE 

CALL PARTA( I »K»LtD) 

23 WH(L»K)*P1(K-1) 

20 WM(L.K)=AM1(K) 

IF(L-l) 10.10.9 

9 DO 3 K=1.N 

DEL=(WM(L.K)-WM(L-1.K) )/WM(L.K) 

IF(MPRNT)10.18.19 

19 WRITE(61.91)K.WM(L.K> .WM(L-l.K) .DEL 

18 IF(ABS( (WM(L.K)-i«/M(L-l.K) )/WM(L.K) ) -5 .OE-03 ) 3 . 3 . 10 

3 CONTINUE 
GO TO 12 

10 IF(L-3)15.16.16 

16 DO 17 K=3.N 

17 P1(K-1)*WH(L.K)-(WM(L.K)-WM(L-1.K) ) * ( WH ( L .K ) -WH ( L“l .K ) )/(WM(L.K) 
12.0*WiM(L-l.K)+WM(L-2.K) ) 

DO 24 K=1.N 

24 AMl(K)*(WM(L.K»+WM(L-l.lO )/2.0 
IF (MPRNT) 15.15.22 

22 WRITE(61.90)L.L.AM1(1) .Aivil(2) .AMK3) .AMI (4) .P1(2).P1(3» 

15 CONTINUE 
GO TO 12 

11 EM=0.0 
EMX=0.0 

DO 14 IP=1.N 
CMK IP)=+0.0 
14 ALl(IP)=+0.0 

12 CALL OUTPUT( I .L) 

WRITE(61.555) 

13 CONTINUE 
STOP 

555 FORMAT (IHO) 

90 FORMAT(2I4.8E15.8) 

91 FORMAT( I4.5E15.8) 

END 
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SUBROUTINE PARTB ( K1 ♦ LI » I) 

DIMENSION PM(41) 

COMMON MPRNT 

COMMON NH»N»NT*NP»GC»OM»TC»PCtG»TT( AO) »HT(AO) »FMT(AO) fFMF(AO) »W»NO 
COMMON CPF( A0*5 ) » A ( S ) »HK ( 5 ) »PL ( D * VO ( 5 ) t AL ( S ) »CCF ( 5 ) » PAl ( S ) »PZ( 1) 
COMMON PUS) »R1(5) »AM115> ♦Tl.CMUS) »ALl(b) 

COMMON P2(5) .R2<5) tAMaiS) tT2»CM2(lj) »AL2(b) 

COMMON WH(A0»5) »TP(A0) *OP(AO) . EP . EK » EM ♦ EMX . PA » T A » KA »C A » ALT » FM 

COMMON WM(AO»5) »OT(AO) »U 

K=K1-1 

NU=0 

NS = 0 

DO 1 L*1»A0 
IF(MPRNT-1)2»2»10 

10 WRITE(61»555) 

WRITE(61.100)L»K 

2 IF(L-2)6»3»5 

6 PM(L)=P1(K1) 

GO TO 5 

3 PM(L)=.998#PM(L-1) 

5 P1(K1)=PM(L) 

TO»Pl(ia)/(Rl(Kl)*GC*9.8) 

CALL ENTR(K»SV) 

PM(L)=P1(K1 ) 

AMI (K1 )=AM1 (K.) 

HPO=EM 

TL=TO/(l.+(G-l.)*EM *EM /2.) 

CALL FRICTIAMKK ) » A ( K ) t TL . FK ♦ V I S # RN ) 

FF=G*(FR*PL(K)/(2.*HR(K) )+W) 

CALL PIPE(FF»CE»SV»HP»TO.TX) 

IF(MPRNT-l) 11»11»12 

12 WRITE(61»100)L»K»L1»P1(K) . EP * EM. EMX *PA1 ( K) » AMI ( K ) . AMI (K 1 ) .HP 
WRITE(61»100)L.K1.L1»P1 (K1 ) »R1 (K1 ) »P2(K1 ) »P2(<) 

11 DE=P1(K)-PA1(K) 

DEM=HPO-HP 

CALL ITER(NU.NS»PM(L) »DE ♦ DEM » PM ( L+1 ) .MSB.MSO) 

IF(MPRNT-l) 13.13.14 

14 WRITE(61.101)L.K.AM1(K1) . PM ( L ) . PM { L+ 1 ) .DE.DEM 

13 IF(EMX-1,)17.7.7 

17 IF(MSB)4.8.4 

4 IF(ABS(DE)-(1.0£-07*P1(K1) ) )8,8.1 

7 IF(DE+(2.0E-07*P1(K1) ) )1.9.9 

9 IF(ABS(DEM)-1.0E-04*EM)8.8.18 

18 IF (MSO) 1.8.1 
1 CONTINUE 

8 WM(L1.K)=AM1(K) 

WM(L1.K.1)=AM1(K.1) 
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CM1(K1)*EM 
CM1(K)=EMX 
IF(MPRNT) 15»15»16 

16 MRITE(61»102)LltL(MSb»MSUiAMKK) *P1(K1) fPZCKI) 
15 RETURN 

101 FORMAT(2I4.6E15.8) 

100 FORMAT(3I4»8E15.8) 

102 FORMAT(4I4,8E15.8) 

555 FORMAT(IHO) 

END 


4 
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SUBROUTINE INPUT 
COMMON MPRNT 

COMMON NH»N»NTtNP#GC»OM»TC»PCfG»TT(40) fHT(40) *FMT(40) t FMF ( 40 ) # W tNO 
COMMON CPF « 40 1 5 ) »A( 5) »HR( 5) »PL ( 1 ) t VO ( 5 ) » AL ( 5 ) # CCF ( 5 ) ♦ PAl ( 5 ) *PZ( 1 ) 
COMMON P1(5)»R1(5) .AMK5) »T1»CM1(5) iALl(5) 

COMMON P2(5) »R2(5) »AM2(5) »T2»CM2(5) .AL2(5) 

COMMON WH(40.5) ♦TP(40» fOP(40) . £P t ER t EM # EMX t PA » TA »RA tCA » ALT » FM 
COMMON WM(40.5) »QT(40) »Q» AM3 ( 5 ) t P3 ( 5 ) 

READ(60»100)NHtNP»MPRNT 

READ(60»101)OM»TC»PC.GC»G 

VC=7.7*SQRT (OM)*(TC**(“l./6« ) ) * I PC** ( 2 • / 3. ) )*1.0E-07 

WRITE(61*126)GC»OM.TC.PC 

WRITE(61»128)VC 

DO 1 1=1 »Nh 

1 READ(60»101) TT( I ) »HT( I ) »FMT( I » »QT( I ) 

DO 2 J=1»NP 

2 READ(60»101)FMF( J) »CPF ( J . 1 ) #CPF J J. 2 » tCPF ( J »3 ) »CPF ( J t4 ) 

READ(60f 100)NO»NT»N 

WRITE(61tl30) 

READ(60fl01) A( 1 ) fHR( 1) »PL( 1 ) tRK 1) tPH 1 ) 

WRITE(61»124)A( 1) »HR{ 1) »PL{1) 

WRITE(61»108) 

DO 5 I=2»N 

READ (60 *101) A( I ) »VO( I ) »R1 ( I ) *P1 ( I ) *AM1( I ) 

READ (60 *101) AL( I) *HR( I ) »CCF( I ) »AL1( I ) «AM2( I ) 

READ(60»101) P2(I) 

R2( I )=R1( I ) 

P3( I )=P2( I ) 

AM3( I )=AM2( I ) 

P2( 1 )=P1( I) 

AL2( I )=AL1( I ) 

AM2( I )=AM1( I ) 

WRlTE(61*123)It A(l)»VO(l)»Rl(I)»Pim 
5 WRITE(61*127) AL ( I ) »HR ( I ) .CCF ( I ) 

AM3(1)=AM2(2) 

AM2(1)=AM1(2) 

EP=P1(2) 

ER=R1(2) 

READ(60*101)W 

WRITE(61»108) 

RETURN 

124 FORMATdOH PIPE AREA. E15 . 8 . lOH HYOR . RAO. .E15 . 8 . 8H PIPE L..E15.8) 
130 FORMAT! 13H PIPE SECTION) 

101 FORMAT! 6E15. 8 ) 

100 F0RMAT(9I3) 

128 FORMATdSH CR I T . VI SCOS I T Y .E13 .6 .6H KG/MS) 

126 FORMATdlH GAS CONST . * E13 .6 .3H M/0K.7H MOL W..E13.6.9H CR. TEMP. .El 
13.6.3H 0K.9H CR.PRES. *E13 .6 *4H ATM) 

127 FORMAT(22H LEAK AREA.E15.8.3H M2.10H HYDR.RA0.E13.6.2H 

1 M.16H RAD.CLEAR/FL.L. .E13.6) 

125 F0RMAT(9H COMP. NR. . I 3 . lOH OR I F . AREA . E 15 . 8 . 3H M2.5H VOL. .E15.8.3H M 
13.8H DENSITY. E15. 8. 8H KGS2/M4.7H PRES. E15.8.6H KG/M2 ) 

108 FORMAT (IHO) 

END 
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SUBROUTINE TRAJ 
PA»AMBIENT PRESSURE 
TA=TEMPERATURE 
RA*DENSITY 
CA»SPEEO OF SOUND 

PAl (N)«PRESSURE OUTSIDE PIPE AND LEAKS AT TIME T1 
DIMENSION WKAO) 

COMMON MPRNT 

COMMON NH*N*NTiNP*GC»OM*TC*PCfG*TT(AO) »HT(A0) »FMT(A0) »FMF( AO ) ♦WtNO 
COMMON CPF { AO ♦ 5 > ♦ A < 5 > *HR ( 5 ) »PH 1 > »V015) »AL(5) #CCF(5) tPAK 5) »PZ( 1 ) 
COMMON P1(5)»R1(5)*AM1(5)»T1*CM1(5) *AL1( 5) 

COMMON P2( 5) *R2 ( 5) »AM2< 5) iT2*CM2 <5) iAL2( 5) 

COMMON WH(A0»5) »TPIA0) *OP(AO) * EP .ER »EM » EMX *PA »TA»RA tCA» ALT » FM 

COMMON WM(A0*5) *UT(A0) .0 

CALL INTER! 1*NH* 1.1. TTfHT»Tl* ALT) 

CALL INTER! l.NH. 1.1. TT.FMT.Tl.FM) 

CALL INTER! l.NH. 1.1. TT.OT.Ti.Q) 

DO 1 I-l.N 
DO 2 J»1.NP 

2 Wl! J)=CPF! J.I ) 

CALL INTER! 1 .NP . 1 . 1 . FMF . W1 .FM. PAl ! I ) ) 

1 CONTINUE 

HG*ALT*6.378178E+06/ !6.378178E+06+ALT ) 

CALL ATMOS!HG.PA.TA.RA.CA) 

DO 3 1*1. N 

3 PAl ! I )«PA + RA*PA1 ! I )*! ! FM#CA ) **2 ) /2 . 0 
RETURN 

END 
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SUBROUTINE PI PE ( FR »CX »CE »H » TO ♦ TX ) 

C SUBROUTINE FOR CALCULATING PIPE FLOW WITH CONST. HEAT ADDITION 
C Q = TOTAL HEAT INPUT 
C PL » PIPE LENGTH 

C TO = TOTAL TEMPERATURE AT ENTRANCE OF PIPE 

C FR = FRICTION TERM (4F)/D 

C G = SPEC. HEAT RATIO 

C GC = GAS CONSTANT 

C EN = ENTRANCE MACH NUMBER 

C PE * PRESSURE AT ENTRANCE OF PIPE 

C EX = EXIT MACH NUMBER 

C PX * PRESSURE AT EXIT 

C RE » DENSITY AT ENTRANCE 

C RX = DENSITY AT EXIT 

C CE SPEED OF SOUND AT ENTRANCE 

C CX SPEED OF SOUND AT EXIT 

C H = ENTRANCE MACH NUMBER FOR CHOKED FLOW. EMX»1 
C PL * PIPE LENGTH 

DIMENSION T(21) »A(21) »FX(51) »FM(51) »TM(51) 

COMMON MPRNT 

COMMON NH»NZ.NT.NP»GC»0M,TC»PC»G»TT(40) .HT ( 40 ) .FT ( AO ) »FMF( 40 ) .W.NO 
COMMON CF(40»5) »AA(5) »HR( 5) »PL( 1 ) »VO(5) .AL(5) »CCF( 5 ) »PA1 ( 5 ) »PZ( i ) 
COMMON Pl( 5) »R1 (5) .AMKS) »TI »CM1(5) .ALKS) 

COMMON P2(5).R2(5) .AM2(S) »T2»CM2(lj) »AL2(5) 

COMMON WH(40»5) »TP(40) »OP(40) .EP.ER.EM.EMX.PA.TA.RA.CA.ALT.FMM 

COMMON WM(40»5) .QT(40) »Q 

N»6 

M*20 

IJ=M+1 

A1=M*N 

DX=1./A1 

CP*G*GC/(G-1. ) 

DT=Q/CP 

JE*N+1 

B=EM*EM 

Gl=(G-l.)/2. 

G2=G-1. 

G3=(G+1. )/2. 

T( 1)=T0 
TM(1)=T0 
A( 1 )=EM**2 
FX( 1)=0. 

FM( 1)=A( 1) 

CALL PIPEA( JE.M.N. ME .T.TM.A. FX.FM.DX.DT.SLP.AM.T01.fr) 

TERM=ME 

X0=TERM/20. 

TX=T(ME+1) 

EMX=FM(ME+1 ) 

EX=EMX 
MA = 0 
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IF(MPRNT-2)3.3t5 

5 WRITE(61»90) JE»M»N»ME»MA»SLP#X0.TX»EMX»EX#T01tDT 

90 F0RMAT(5I3»8E15.8) 

3 IF(SLP-.0020»65*65»66 
66 TERM=ME-1 

X0*TERM/20. 

AMaA(ME) 

T01*T(ME) 

TM(ME>*T(ME) 

65 IF(XO-l. )7.9.9 

7 N*20 
A1*N 

DM2*( l.-AM) /20. 

DM1=DM2/A1 
A( 1 )=AM 
T( 1 )=T01 

CALL PIPEB( JEtM»N»MA»ME»T.TM»A»FX»FM»0T»0Ml»DM2»TX»FR»X0) 
9 CONTINUE 
M1=ME+MA 
KK»1 

TM(M1)=TX 
IF(EMX-1. ) 1»2»2 

1 FX(M1)»1.0 
FM(M1)»EMX 
EXA»SQRT(EMX» 

H=EM*SQRT(1«/EXA) 

GO TO 4 

2 EMX=1.0 
EXAal.O 

H=EM*SQRT(FX(M1) ) 

4 DST=TX*( 1.+G1*EM*EM)/(T0*(1.+G1*EMX) ) 

PI ( 1 ) =EP*EM*SQRT ( OST ) /EXA 
R1(1)=ER*P1(1)/(DST*EP) 

CX = SQRT(G*P1( 1) /Rl( U ) 

EMX=SQRT(EMX) 

IF(MPRNT-2)6»6»8 

8 WRITE(61»91)M1»KK»EMX.FX(M1) fTMIMl) »EXA»H»P1{1) »R1(1) 

91 FORMAT(2I3»8E15.8) 

WRITE(61»91)M1»KK»TX»T0»DT*PA1 ( 1 ) 

6 RETURN 
END 
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SUBROUTINE P I PEA ( JE »M i N *ME » T » TM» Af FX »FM i DX »UT f SIP * AM» TUI iFR ) 
DIMENSION T(21)*A(21)»FX(51)iFM(51)»TM(51) 

DIMENSION Y(91)»S(91)»Tl(91)iSl(91)»BM(91) 

COMMON MPRNT 

COMMON NH»NZ.NT»NP»GC»OMiTC»PCiGfTT(AO) »HT ( AO ) »FT ( AU ) »FMF(A0) »W»NO 
COMMON CF( AO*5) »AA(5> *HR( 5) »PL{ 1 ) *VO(5) iAL( 5 ) *CCF( 5) i PAl ( 5) *PZ( 1 > 
COMMON Pl(5)»Rl(5)*AMl(b)tTI iCMl (5) »AL1 ( 5) 

COMMON P2(5) .R2(5) *AM2(5) .T2 #CM2 ( 5 ) » AL2 ( 5 ) 

COMMON i*/H(AO»5) *TP(A0) *OP(AO) i EP i ER * EM i EMX »PA * TA . KA i CA * ALT ♦ FMM 
COMMON WM(AOi5) »QT(AO) *Q 
TOL*1.0E-06 
I J*M+1 

Gl* (G-1. ) n* 

G2»G-1. 

G3*(G+1. )/2. 

DO 1 I*1»M 

A8 A1«A( I )*( l.+G*A( I) )*{ l.+Gl*A( I ) )»DT/ ( T ( I ) * ( 1 .-A (I ) ) ) 

A2*G*( A( I )#*2 )*( l.+Gl#A( I ) )*FR/ ( 1.0* ( l.-A( I ) ) ) 

DO A2 J=1»JE 
X»J-1 

Y( J) = (A1+A2)*X*DX+A( I ) 

S( J)=Y( J) 

A2 Tl( J)»T( I >+DT*DX#X 
T( I+1)=T1( JE) 

TM( I+l )=T( I+l ) 

A6 DO 2 J-l.JE 

A1=S( J) *( l.+G*S( J) )*( l.+Gl*S(J) )*DT/ ( 1.0*T1 ( J )*( l.-S( J) ) ) 
A2=G*(S(J)»*2)*(1.+G1*S(J) )*FR/(1.0*(1.-S( J) ) ) 

2 Sl( J)»(A1+A2) 

CALL INTI JE»DX»S1*BM) 

DO 3 J=1»JE 
3 S< J)=A( I )+BM( J) 

IF(S( JE) )70»71.71 

71 IF(S( JE)-.95)72.72.70 
70 NL=JE+1 

GO TO A9 

72 DO 43 J=1»JE 
IF(MPRNT-5)A.A»5 

5 SSJ=SQRT(S( J) ) 

SYJ=SQRT(Y(J) ) 

WRITE(61.90)I »J. JE.SSJ.SYJ.TOL 
A IF(ABSF(SQRTF(S( J) )-SQRTF( Y( J) ) ) -TOL ) A3 . A3 . AA 
A3 CONTINUE 
NL=JE+1 
GO TO 20 
AA NL=J 

20 DO 35 J=1.JE 
35 Y(J)=S(J) 

SLP=S1 ( JE)*DX 
IF( I-I J) 36.3A.3A 
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36 DO 45 J»1#JE 

IFiSl ( J)*DX-.0020»45t45#33 
45 CONTINUE 
GO TO 34 

33 IF(JE-60)51»51*49 

51 N=N+4 

JE*N+1 
A1*M#N 
OX«l,/Al 
GO TO 48 
49 U»I 

34 A(I+1)»S(JE) 

ME*I 

AM=A( I+l) 

T01=T( I+l) 

IF(NL-( JE+1) )46»37t37 
37 A1=I 

FX( I+l)=Al/20. 

FM( I+1)=A(I+1) 

TM( I+1)=T( I+l) 

IF( I-IJ) 1»54»54 
1 CONTINUE 
54 RETURN 

90 FORMAT(3I4»8E15.8) 

END 


Jt 
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SUBROUTINE P I PEB UE »M .N tMA»ME t T » TM. A . FX * FM . DT * DM1 . DM2 » TX »FR »XO ) 
DIMENSION Y(91)»Tl(91)»Z(91)tS(91)iSl(91)*BM(91)*R(91)iP(51)»C(51) 
DIMENSION T(21) »A(21) »FX(51) *FM ( 51 ) » TM ( 5 1 ) 

COMMON MPRNT 

COMMON NH»NZ»NT»NP*GC»OM.TCiPC*G»TT ( AO ) * HT ( AO ) » FT ( AO ) »FMF( AO) tW»NO 
COMMON CF( AO»5) *AA( 5) »HR( 5 ) fPL( 1 ) »VO(5) »AL( 5) »CCF( 5 ) »PA1( 5) »PZ( 1 ) 
COMMON Pl( 5) *R1 ( 5) *AM1 ( 5) »TI .CMl (5 I »AL1 1 5) 

COMMON P2 ( 5 ) »R2 ( 5 ) *AM2 ( 5 ) .T2*CM2 ( 5 ) *AL2 (5) 

COMMON WH(A0.5) »TP(AO) .OP(AO) » EP * EK . EM » EMX f PA i T A t K A »C A i ALT f FMM 
COMMON WM(A0.5) »QT<A0) *Q 
TOLsl.OE-06 
Gl= (G-1. ) /2. 

G2*G-1. 

G3=(G+1, ) /2. 

C( 1 )=XO 
DO A I*1»M 

53 Al»( l.-A( I ) )/(A( I )*( l.+Gl»A( I ) ) ) 

A2»l./( ( l.+G*A( I ) )*DT/T( I l+G*A( I )*FR) 

SLP*A1*A2*DM1 

IF(N-60)62»62»61 

62 IF(N-6)60»63»63 

63 IF(A1*A2*DM1-. 0025)39. 38»A0 
39 IF(A1*A2*DM1-.0015)A1»A1»38 
A1 N«N-2 

JE»N+1 

IF{N-6)60»52.52 

60 N*6 
JE»7 

GO TO 38 
AO N=N+A 
JE=N+1 

IF(N-60)52»61.61 

61 N=60 
JE = 61 

GO TO 38 
52 A1=N 

DM1=DM2/A1 
GO TO 53 
38 DO 5 J»1.JE 
X*J-1 

IF(J-l) 18.18.19 

18 Y( J)=C( I ) 

Tl( J)*T( I ) 

Z( J)=Y( J) 

GO TO 5 

19 Y( J)=Y( J-1)+A1*A2*DM1 

Tl( J)=T1( J-1)+DT*(Y( J)-Y( J-1) ) 

Z( J)=Y( J) 

5 S( J)=A( I )+DMl*X 
A( I+l )=S( JE) 
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lU DO 8 J=1»JE 

Al»(l,-S< J) >/(S( J)*(1.+G1*SU) ) ) 
A2=l./( ( l.+G*S( J) )*DT/T1( J)+G*S( J)*FR) 
8 Sl( J)*(A1*A2) 

CALL INT(JE»OMl»Sl»BM) 

DO 10 J*1»JE 
2(J)»C( n+BM( J) 

IF( J-l)16»16.17 

16 Tl( J)»T( I ) 

GO TO 10 

17 Tl( J)=T1( J-1)+DT*(Z( J)-Z( J-1) ) 

10 CONTINUE 

DO 11 J»1»JE 
IF(MPRNT-5) 1»1»2 

2 WRITE (61 *90) I * J ♦ JE »Z ( J ) » Y U ) »TOL 
1 1F(ABSF(Z(J)-Y( J) )-TOL>ll»ll»12 

11 CONTINUE 
NL*JE+1 
GO TO 21 

12 NL=J 

21 DO 13 JaltJE 

13 Y(J)=Z(J) 

C( I+1)«Z( JE) 

T( I+1)»T1( JE) 

IF(NL-( JE+1) ) 1A»6»6 
6 MA=I 
M2=ME+1 
FX(M2)»C( I+l) 

FM(M2)»A( I+l) 

TM(M2)=T( I+l) 

IF(C( I+l)-l.)22»26»26 

26 R(l)»l. 

CALL INTER! 1»JE»1»1»Z»S»R»P) 

EMX*P( 1) 

CALL INTER! 1»JE.1»1»Z»T1.R»P) 

TX=P! 1) 

GO TO 27 

22 IF!A! I + D-l. ) 15»28.28 
28 EMXsl.O 
TX=T! I+l) 

GO TO 27 

15 IF! I-M)4»28»28 
4 CONTINUE 

27 RETURN 

90 FORMAT!3I4»8E15.8) 

END 
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SUBROUTINE I TER ( NU »NS fPiVI .DE »OEM» PMl »MSUB»MSON ) 

DIMENSION DPX(40) »AMA(40) »DPE(40) »AMB(40) 

COMMON MPRNT 

COMMON NH»N*NT»NP.GC*OM»TC»PC»G»TT(40) »HT(40) »FMT(40) »FMF(40) »W»NO 
COMMON CPF(40»5 ) » A ( 5 ) *HR ( 5 ) ♦ PL ( 1 ) ♦ VO ( 5 ) ♦ AL ( 5 ) »CCF ( 5 ) ♦ PAl ( 5 ) »PZ( 1) 
COMMON PI ( 5) »R1 (6) tAMUS) .TltCMl (5) »AL1( 5) 

COMMON P2( 5) tR2 (5) »AM2(5 ) »T2»CM2(5) »AL2( b) 

COMMON WH(40»5) ♦TP(40) »OP(40) »EP*ER»EM»EMX 

MSUB=1 

MSON=l 

IF(EMX-1. )1»3»3 

1 NU*NU+1 
ANU*NU 
OPX(NU)=DE 
AMA(NU)=PM 
IF(NU-2)6»2»2 

2 AA=ABS(AMA{NU)-AMA(NU-1) ) 

IF(AA-1.0E-07)9»12»12 

12 PM1=AMA(NU)-DPX(NU)*(AMA(NU)-AMA(NU-1) ) / ( DPX ( NU ) -DPX ( NU“1 ) ) 

GO TO 5 

9 MSUB=0 

IF( ABS(DE)-( 1.0E-07*P1 (2 ) ) )8»8.15 

15 IF{ABS(DEM)-2.0E-06)16.16»8 

16 MSUB=1 
EMX=1.0 

3 NS=NS+1 
ANS=NS 
DPE(NS)=DEM 
AMB(NS)=PM 
IF(NS-2)8*4»4 

4 AB=ABS(AMB(NS)-AMB(NS-1) ) 

IF(AB-1.0E-07) 13»14.14 

13 MS0N=0 
GO TO 8 

14 PM1=AMB(NS)-DPE(NS)*(AMB{NS)“AMB(NS-1) ) / ( DPE ( NS ) -DPE ( NS“1 ) ) 
IF(PM1-AMB(NS) ) 10»11»10 

11 PM1=AMB(NS)*(ANS+0.010)/ANS 

10 IF(ABS(DEM)-1*0E-04*EM)7»7»5 

7 IF(DE+1.0E-07)1»5.5 

6 PM1=(PM+P2(2) )/2.0 

GO TO 5 

8 PM1=.9995*PM 

5 RETURN 
END 
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SUBROUTINE OUTPUT ( INT »L ) 

COMMON MPRNT 

COMMON NH.N»NT#NP»GC»OM»TC»PC»GiTT(AO) »HT ( 40 » #FMT ( 40 ) #FMF(40) »W»NO 
COMMON CPF(40»5) .A ( 5 ) .HR ( 5 ) »Plt ( 1 ) »V0 ( 5 ) » AL ( 5 ) »CCF( 5 ) »PAl ( 5 ) »PZ( 1) 

COMMON Pl(5) »R1(5) »AMK5) .T1.CMK5) »ALl(5) 

COMMON P2(5).R2(5) .AM2(5) .T2.CM2(5) »AL2(6) 

COMMON WH(40.5> »TP(40 ) »OP(40) . EP »ER » EM. EMX »PA » TA »RA »CA* ALT » FM 
COMMON WM(40»5) »QT(40 ) »Q»AM3( 5 ) .P3( 5 ) 

WRITE(61»50)T1.L 

50 FORMATdlH TIME E12.5.11H NR. ITER. 13) 

WRITE(61.60) ALT.FM.PA.TA.KA 

60 FORMATdlH ALTITUDE E12.5.11H MACH NR. E12.5.11H AMB.PR. 612.5 
l.llH AMB.TEMP E12.5.11H AMB.OENS E12.5) 

WRITE(61»61 )EP#ER»£M»AM1 d) 

61 FORMAT(34H COMPARTM. PIPE ENTR.PR. E12.5.11H ENTR.DENS.El % 

12.5.11H ENTR.M.NR.E12.5.11H ENTR.M. FL. E12. 5 ) 

WRITE(61.62»Pld) »Rld) .EMX.AMKl) 

62 FORMAT(34H EXIT PR. E12.5.11H EXIT. DENS El 

12.5.11H EXIT M.NR.E12.5.11H EXIT M.FL.E12.5) 

WRITE(61.65)Q 

65 FORMAT! 34H HEAT AD. E12.5) 

DO 1 I=2»N 
0PN=P1( I)-PA 

WRITE(61»63)I»Pld).Rld) »CM1 ( I ) »AM1 ( I ) 

1 WRITE(61.64)PA1 ( I ) .ALld ) »DPN 

63 FORMAT(12H COMPARTM. 13.19H PRESSURE E12.5.11H DENSITY 

1 E12.5.11H MACH NR E12.5.11H MASS FLOW E12.5) 

64 FORMAT(34H LEAK PR. E12.5.11H LEAK M.FL.El 

12.5.11H PRES.DIFF.E12.5 ) 

IF( INT-NT)4»6.6 
6 N01=NT+1 

9 NT1=201 

WRITE(62.80)N01.NT1»N 

WRITE(62»81)A(l)»HRd)»PLd).Rld).Pld) 

DO 2 I*2.N 

WRITE(62»81 )Ad ) »VOd ) »R1 d ) .Pld ) .AMld ) 

WRITE(62.81)AL( I ) »HR( I ) .CCFd ) .ALld ) .AM2( I ) 

2 WRITE(62.81)P2( I ) 

WRITE(62.81)W 

GO TO 10 

4 GO TO (8.10) .SSWTCHFd) 

8 N01=INT+1 

GO TO 9 

10 DO 3 I=1.N 

AM3( I )=AM2( I ) ^ 
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P3( I )=P2( 1 ) 

AM2( I )=AM1( I ) 

IF( INT-1)5.7*5 
7 AM2 ( I ) »AM1 ( I ) +1 #0E-06 
5 P2 ( I )*P1 ( I ) 

R2( I )*R1( I ) 

CM2 ( I >=CM1 ( I ) 

3 AL2( I )=AL1( I ) 

AM2 (N+1 ) *0.0 
AM1(N+1)=0.0 
EP*.99*EP 
RETURN 

80 FORMAT(9I3) 

81 FORMAT ( 5E15. 8 ) 

END 
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